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Abstract 

Our goal is to understand primary accretion of the first planetesimals. Some exam- 
ples are seen today in the asteroid belt, providing the parent bodies for the primitive 
meteorites. The primitive meteorite record suggests that sizeable planetesimals formed 
over a period longer than a million years, each of which being composed entirely of an 
unusual, but homogeneous, mixture of mm-size particles. We sketch a scenario that 
might help explain how this occurred, in which primary accretion of 10-lOOkm size 
planetesimals proceeds directly, if sporadically, from aerodynamically-sorted mm-size 
particles (generically "chondrules" ) . These planetesimal sizes are in general agreement 
with the currently observed asteroid mass peak near 100km diameter, which has been 
identified as a "fossil" property of the pre-erosion, pre-depletion population. We extend 
our primary accretion theory to make predictions for outer solar system planetesimals, 
which may also have a preferred size in the 100km diameter range. We estimate for- 
mation rates of planetesimals and explore parameter space to assess the conditions 
needed to match estimates of both asteroid and Kuiper Belt Object (KBO) formation 
rates. For parameters that satisfy observed mass accretion rates of Myr-old protoplan- 
etary nebulae, the scenario is roughly consistent with not only the "fossil" sizes of the 
asteroids, and their estimated production rates, but also with the observed spread in 
formation ages of chondrules in a given chondrite, and with a tolerably small radial dif- 
fusive mixing during this time between formation and accretion. As previously noted, 
the model naturally helps explain the peculiar size distribution of chondrules within 
such objects. The optimum range of parameters, however, represents a higher gas 
density and fractional abundance of solids, and a smaller difference between keplerian 
and pressure-supported orbital velocities, than "canonical" models of the solar nebula. 
We discuss several potential explanations for these differences. The scenario also pro- 
duces 10-lOOkm diameter primary KBOs, and also requires an enhanced abundance 
of solids to match the mass production rate estimates for KBOs (and presumably the 
planetesimal precursors of the ice giants themselves). We discuss the advantages and 
plausibility of the scenario, outstanding issues, and future directions of research. 



1 Introduction 



Primary accretion is the stage of growth in which tiny protoplanetary nebula dust grains 
grow into objects of 10-100 km size, such as most asteroids, Kuiper Belt Objects (KBOs), 
and comets. The most well known, traditional approaches to modeling primary accretion 
are incremental growth by simple sticking to ever larger sizes (Weidenschilling 1997, 2000; 
DuUemond and Dominik 2004, 2005; reviewed by Dominik et al 2007; most recently Brauer 
et al 2008) and midplane instabilities of various types, going back to Goldreich and Ward 
(1973) {of. reviews by Cuzzi and Weidenschilling 2006; henceforth CW06, and more recently 
Chiang and Youdin 2009). 

Several important clues as to the nature of primary accretion, which can help us assess 
these different hypotheses, are to be found in primitive meteorites and asteroids (discussed 
in more detail in section 2.1). The most primitive chondritic meteorites display a character- 
istic texture: predominance of mm-sized, once-molten silicate chondrules, metal grains, and 
refractory oxide particles, each surrounded by fine-grained dust rims and all embedded in a 
granular matrix. The size distribution of the chondrules in all classes of chondrite is quite 
narrow and nearly universal in shape, but with a mean size distinctive of each class. At 
least two entire chondrite classes are each thought to derive from only one or two planetesi- 
mals, roughly 100 km in size and originally composed largely of chondrules with very similar 
properties. This ubiquitous and unusual texture is surely telling us something important 
about primary accretion, but there is no explanation for it at present. The Myr duration 
of meteorite parent body formation as revealed in isotopic age-dating, and the prevalence of 
unmelted asteroids, suggest that primary accretion went on for a long time (section 2.1). 

The observations suggest that primary accretion was inefficient, and took a long time 
to complete (CW06; Cuzzi et al 2008, henceforth CHS08; see section 2.1). If the nebula 
were nonturbulent, as required for traditional midplane instabilities to play a role, particles 
settle into a dense midplane layer and growth by incremental accretion comes to completion 
too quickly - numerous 100km planetesimals and even lunar-size objects form in 10^ years 
(Weidenschilling 2000), all of which would melt due to short-lived radionuclides such as ^^Al. 
Nebula gas turbulence can frustrate primary accretion if simple "incremental accretion" 
stalls at roughly dm to m-size in turbulence, depending on gas density (the so-called "m- 
size barrier"; see Cuzzi and Weidenschilling 2006, Dominik et al 2007, Brauer et al 2008). 
Moreover, recent work has identified a new "km-size barrier" for incremental accretion in 
turbulence (Ida et al 2009). A challenge for primary accretion in turbulence is to leapfrog 
not only the meter-size barrier, but perhaps also the km-size barrier - and create 10- 100km 
asteroids entirely from "chondrules" with similar properties. It is in this sense that the first 
planetesimals might indeed have been 10-lOOkm in diameter. If this happens in a temporally 
extended fashion, nebula chemical and physical properties can change slowly, perhaps helping 
explain the variable chemical and isotopic properties of chondrites {eg Cuzzi et al 2005). 

In previous work we have emphasized intriguing connections between these properties of 
primitive meteorites and asteroids, and the general scenario we present here. We have shown 
how well-sorted, chondrule-sized mineral particles are concentrated, by orders of magnitude, 
into dense zones in weak nebula turbulence (sections 2.2-2.3). This turbulent concentration 
can explain the characteristic size and size distribution of chondrules in a natural way. We 
developed a cascade model of the statistics of dense zones and their correlation with gas 
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vorticity, which incorporates the effects of particle mass loading on the gas and predicts 
the fractional volume of particle-rich zones which can evolve directly into objects with some 
physical cohesiveness. Here we derive threshold conditions (combinations of the density and 
lengthscale of particle clumps, and the density, pressure gradient, and local vorticity of the 
gas) which allow dense clumps to proceed to become actual planetesimals. Combination of 
these thresholds with our cascade models leads to a prediction of the relative abundance of 
primary planetesimals as a function of mass - their initial mass functions - and even (with 
uncertainties) their production rate (sections 3.1-3.4). 

In this paper we explore, in a preliminary way, primary accretion initial mass functions 
(IMFs) at two disparate locations in the early solar system. These predictions may be 
compared with both the known asteroid size distribution (which is thought to be a "fossil" 
representing the actual size of primary planetesimals), and also with (limited) knowledge 
for Kuiper Belt Objects in the 30AU region. Under different assumptions regarding nebula 
properties, we estimate not only the characteristic planetesimal size or mass which results, 
but also the planetesimal formation rate, which can itself be compared with crude estimates 
in the asteroid and KBO regions (section 3.4). Our IMFs are consistent with previous sug- 
gestions that "asteroids were born big" (Bottke et al 2005; Ida et al 2008; Morbidelli et al 
2009a; WeidenschiUing 2009). The scenario we envision for primary accretion, based on tur- 
bulent concentration, might occur continuously - and inefficiently - over an extended time, 
but when it does occur it is highly selective as to constituents and bypasses the problematic 
meter-size range (and even the km-size range) entirely, leading directly to 10-lOOkm size 
objects composed of aerodynamically sorted particles. An independent study along these 
lines has also been done by Chambers (2010). A different scenario has been advanced to 
explain direct growth to 100km or larger diameter bodies, starting with meter-size bodies 
(Johansen et al 2007). This alternate pathway occurs in environments similar to that de- 
scribed here, and could proceed simultaneously (see sections 2.2 and 4 for more discussion). 
It will become apparent that current uncertainties in both the observations and the theory 
render our scenario more of a suggestive roadmap for, rather than an exhaustive explanation 
of, primary accretion. 

2 Background 

2.1 Clues from meteorites, asteroids, and KBOs: 

Meteorites: Several different isotope systems (Al-Mg and Pb-Pb primarily) testify that the 
bulk of chondrites (their dominant iron-magnesium-silicate chondrules and matrix) was last 
processed in the nebula 1-3 Myr after the formation of the oldest, highest-temperature 
minerals (the refractory Calcium- Aluminum-rich Inclusions or CAIs) found in the same me- 
teorites (Russell et al 2006; Kita et al 2000, 2005). The rare, even later-forming CH and 
CB chondrites probably resulted from an entirely different process, in an entirely differ- 
ent environment (Wasson and Kallemeyn 1990, Krot et al 2005). Yet, some parent bodies 
apparently accreted and melted nearly contemporaneously with CAIs, forming achondrites 
and metal cores (Kleine et al 2005, Markowski et al 2007). Primary accretion thus lasted 
several million years, suggesting that it was inefficient. Moreover, isotopic age-dating has 
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recently progressed in accuracy and quantity to the point where several different groups 
find, for several different chondrite classes (carbonaceous and ordinary), that the formation 
ages of chondrules within a given chondrite range over almost 1 Myr (Kita et al 2000, 2005; 
Mostefaoui et al 2002, Sugiura and Krot 2007, Kurahashi et al 2008, Villcnueve et al 2009). 
The nominal two-sigma error bars on individual chondrite ages in the best of these data are 
roughly 0.3-0.4 Myr, so a range of perhaps a half-million years can't be ruled out, but it 
appears from taking these results at face value that a range as short as 10^ — 10^ years is 
unlikely in spite of qualitative thinking in the past that chondrules had to be accreted into 
chondrites "rapidly" after their formation. Cautionary notes have been raised that some or 
all of these apparent age spreads might be the result of mineral-specific parent body alter- 
ation processes {eg. Alexander 2005, section 7.2); it is of vital importance to continue to 
make measurements of this type while addressing questions of alteration because, as we will 
show, they provide powerful constraints on models of primary accretion. 

The texture of primitive chondrites is unusual, and suggests a role for aerodynamical 
effects in most cases (for reviews see Scott and Krot 2005 or Brearley and Jones 1998; 
Cuzzi 2004 and CW06 present more discussion of the evidence for aerodynamical effects). 
The sizes of sihcate and metal particles in the young CH and CB chondrites are counter- 
indicative of aerodynamic sorting, showing the prevalent evidence from normal chondrites 
to be non-trivial. The most primitive chondrites - especially those containing unbrecciated 
"primary texture" (Metzler et al 1992, Brearley 1993) - look like collections of dust-rimmed 
chondrules and other mm-size particles, directly accumulated and merely compressed and 
compacted. Individual constituents of chondrites (chondrules in particular) have a size dis- 
tribution that, while centered at different sizes from class to class, has a not-quite-lognormal 
shape that appears universal (CHPDOl, Teitler et al 2009; sec section 2.2). The H-type 
ordinary chondrite class is believed to derive from a single 80-100 km radius parent body, 
initially composed of a homogeneous collection of similarly well-defined chondrules which 
experienced post-accretional heating, metamorphism, and cooling to different degrees at dif- 
ferent depths (Trieloff et al 2003, Grimm et al 2005). A similar (but less clear) story can be 
told for the L and LL-type ordinary chondrites (Marti and Graf 1992). It's reasonable to 
suspect that chondrite parent bodies may all be large (~100 km) objects, each initially com- 
posed primarily and homogeneously of chondrules (and other associated mineral particles) 
with average chemical, physical, and isotopic properties which are well-defined in any parent 
body, but differ dramatically from one parent body to another. Thus, primary accretion 
may be inefficient, but when it operates, it is highly selective. We return to an assessment 
of the situation in our concluding remarks. 

Asteroids: Most of the S-type asteroids are probably related to ordinary (unmelted) 
chondrites (Binzel et al 2002, Clark et al 2002). This is not to say their interiors were 
never heated, or even partially melted (Elkins-Tanton and Weiss 2009) but there are only 
a few asteroid surfaces manifesting widespread and complete melting, as on Vesta. For 
instance, Sunshine et al (2004) show that in addition to Vesta and the unrelated, but similarly 
differentiated basaltic object 1489Magnaya, three other S-type family parents (17Thetis, 
847Agnia, and 808Merxia) have igneous surfaces. Others of this type might yet be found. 
However, all objects larger than 50 km radius would melt extensively if they accreted earlier 
than 1.5-2.5 Myr after CAIs, because of radiogenic heating by live ^^Al (LaTourrcttc and 
Wasserburg 1998, Woolum and Cassen 1999, McSween et al 2002, Hevey and Sanders 2006). 
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The combination of few thoroughly melted asteroids and many unmelted ones, like the spread 
in meteorite age dates, points to a temporally extended primary accretion process. 

The observed asteroid population shows a distinct mode in the distribution of mass as 
a function of size Figure 1. Bottkc et al (2005) locate the peak of the observed asteroid 
mass distribution at 100km diameter, using cumulative distributions. They make a case 
that this mass peak is not explainable by erosive processes, and instead testifies to an initial 
mass function deficient in smaller objects. Figure 1 shows a differential presentation of 
the asteroid data that suggests the mass peak may he at 140km diameter. For diameters 
larger than 350km, there are two or fewer asteroids per bin, so the details of the distribution 
are highly uncertain; for comparison there are about 50 asteroids in the mass bin at 140 
km diameter. Nevertheless, it is a fact that the asteroid belt mass is dominated by the 
few largest asteroids. This is generally taken as evidence for runaway accretion into even 
larger objects, of which more than 99% have been subsequently removed by size-independent 
dynamical depletion processes (Chambers 2004). The peristence of the 140km bump testifies 
to the vast number of asteroids of these sizes in the prc-dcplction population; the diameter 
at the peak is thought to be an unbiased estimate of the primitive asteroid mass distribution 
at the time dynamical stirring and removal began (presumably at the time the nebula gas 
was removed and/or Jupiter formed; see Bottke et al 2005 or Morbidelli et al 2009a for a 
discussion). Whether the actual primary bodies needed to be just the same size as the current 
fossil population (Morbidelli et al 2009a), or a factor of 3-10 smaller in diameter, incurring 
some subsequent growth before the start of the erosive regime (Weidenschilling 2009) remains 
a subject of debate. Either way, our predictions of the IMF and other physical properties of 
primary bodies provide initial conditions for, and are testable by, models such as these. 

Kuiper Belt Objects: The KBO size distribution, and indeed the entire KBO formation 
scenario, is less well constrained. It is generally agreed that there is a KBO "mass bump" as 
in the asteroid case (figure 1), but the modal peak may lie anywhere between 20 and 100km 
diameter based on the same (magnitude) data, given uncertainties in the observations and 
assumed albedos. Values close to the low end of this range might be ascribed to coUisional 
erosion, for weak objects (Kenyon et al 2008), but values at the high end of this range would 
probably represent "fossil" signatures of the primary accretion process, as in the asteroid 
case (Bottke et al 2005). KBOs come in several dynamical classes, which have different 
size distributions at sizes larger than the mass bump (Bernstein et al 2004, Morbidelli et al 
2009b); of these the most abundant are the "hot" and "cold" classical objects, named for 
their relatively high and low eccentricities and inclinations, respectively. 

It is not even known for certain where the currently observed KBOs were formed. Tra- 
ditional incremental accretion models form them in place (Stern and Colwell 1997, Kenyon 
and Luu 1998). This scenario requires a massive local source population of solids, of which 
more than 99% must be subsequently removed by erosive collisions or by dynamics. Cleanup 
by dynamics alone is slow, unless augmented by local embryo-size objects which are, as in 
the asteroid belt, subsequently lost themselves (Chiang et al 2007). Ford and Chiang (2007) 
explored the excitation of KBOs by local icy embryos which were subsequently lost, with 
mixed results, but detailed studies of dynamical clearing per se in this scenario have not yet 
been done. Kenyon et al (2008) beheve that cleanup by erosion (with planetesimals stirred 
only by Neptune) can remove more than 90% of the bodies by grinding and drag loss of small 
particles (assuming a size distribution with plentiful 1-lOkm bodies); however, removal by 
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Figure 1: A histogram indicating where the bulk of mass Ues in the current asteroids. Diameter bin 
centers are spaced by a factor of two in mass (ie, in D^), starting at the largest asteroid (Ceres) and ranging 
downwards. Diameter boundaries between bins are taken midway between diameters at bin centers. Results 
are shown from several asteroid databases: the IRAS albedo-diameter data posted on the PDS Small Bodies 
node (IRAS-A-FPA-3-RDR-IMPS-V6-0; red squares), a tabulation by FarineUa and Davis (1992, ascribed to 
Tedesco 1989; blue diamonds), and a tabulation by Jedicke et al (2003; green triangles). Poisson statistics 
error bars are indicated. In this representation, most of the mass seems to lie at around 140km diameter. It is 
the conclusion of Bottke et al (2005) that the depletion shortward of 140km diameter cannot be accomplished 
by erosion during post-accretional collisional evolution, but must be a primordial signature. 



erosion would be much less efficient if planetesimals are typically "born big" as implied by 
the larger end of the diameter interpretations of the observations (Bernstein et al 2004). 

The outward dynamical evolution of giant planets by interactions with planetesimals 
(Malhotra 1995, Fernandez and Ip 1996), more recently refined into the so-called "Nice 
model" , has several implications. One is that a massive indigenous population of planetes- 
imals outside of 30AU would lead to greater migration of Neptune than observed, and its 
existence has been questioned on these grounds (Gomes et al 2004). Another suggestion is 
that the current crop of KBOs (of all dynamical types) may have been formed at smaller 
distances - specifically between 16-30AU, and emplaced into their current locations by dy- 
namical processes (Gomes 2003, Levison et al 2008). Supportive connections have been made 
between this emplacement and other related primitive body populations (D-type and Tro- 
jan asteroids) emplaced at the same time and in the same general way (Levison et al 2009, 
Morbidelli et al 2009). In this scenario, there is no "cleanup" problem - the mass emplaced 
into the current Kuiper belt is only about 0.1% of the mass in its 16-30 AU source region 
(and most nebula models contain more than enough total mass in this region). However, the 
observed "cold classical" KBO population is less eccentric than the model results predict; 
also, relative to the hot population, it is four times richer in binaries which might be easily 
disrupted during extended dynamical evolution from closer to the sun (Stephens and Noll 
2006). 

Other issues regard timescales. Some in situ, incremental growth, massive source scenar- 
ios take 10-30Myr to build 10-lOOkm radius KBOs (Kenyon 2002, Kenyon et al 2008), which 
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probably precludes any thermal alteration by short-lived radioisotopes such as ^®A1. On the 
other hand, WeidenschiUing (1997; his figure 12) includes gas-drag augmentation of growth, 
and grows 10-lOOkm objects at 30 AU in <1 Myr; this short accretion timescale would seem 
to predict widespread melting. 

McKinnon et al (2008) note that the emerging diversity of KBO albedos and densities 
(including the apparent differentiation of Pluto and Haumea, formerly 2003EL61, and now 
probably Quaoar as well (Eraser and Brown, 2009) might be hard to explain without short- 
lived isotopes; yet the need to preserve amorphous ice and supervolatiles like CO in other 
objects seems to preclude them. This paradox is reminiscent of the meteoritical arguments 
for drawn-out accretion over a period spanning a little longer than the lifetime of the short- 
lived isotopes. 

Finally, none of the current KBO incremental growth scenarios involve nebula turbulence, 
which slows growth beyond a meter or frustrates it entirely (stalling probably occurs at even 
smaller sizes at these distances, as noted below), and yet it is generally agreed that, if any 
part of the nebula were robustly turbulent, the >20 AU region would be (sec next subsection). 
Even if the meter- (or smaller) size barrier could be overcome in a turbulent environment at 
30 AU, Ida et al (2008) have shown that expected levels of turbulence may excite random 
velocities that render the entire region erosive for 1-lOkm sized objects. No detailed models 
have been run for growth of planctesimals in the 16-30 AU source region, with or without 
turbulence, but at least without turbulence, timescales would probably be faster than in the 
traditional in situ models (WeidenschiUing 1997, Kenyon 2002, Kenyon et al 2008) because 
of higher sohds densities and shorter timescales. 

As in the asteroid belt region, it is plausible that a primary accretion scenario in which 
large planctesimals are created sporadically, over this period of time, and perhaps with a 
different efficiency than incremental growth models, might help resolve some of these KBO 
puzzles. In this paper we will make representative calculations at 30 AU, assuming the 
primordial KBO mass needed between 16-30 AU was about AOM^. Future refinements of 
this preliminary study are discussed in section 3.4.2. 

2.2 Turbulence and particle- gas interactions 

While the ultimate cause and intensity of nebula turbulence remain subjects of debate on 
theoretical grounds (Fleming and Stone 2003; Johnson and Gammie 2005; Turner et al 2007), 
observational arguments suggest it was indeed present at interesting levels throughout the 
primary accretion stage (Dullemond and Dominik 2004, 2005, Dominik et al 2007). The 
most generally accepted (although perhaps not the only) way to drive nebula turbulence 
is the magnetorotational instability (MRI), in which the turbulent intensity is considerably 
higher in the dilute gas of the outer (and upper) nebula than in the terrestrial planet region 
(Turner and Sano 2008). In contrast to the original idea of a "dead zone" near the nebula 
midplanc (Gammie 1996), Turner and Sano (2008) dub the midplane region the "undcad 
zone" because it can be excited in as-yet poorly understood ways by strong turbulence in 
the rarified layers at high altitudes. Moreover, even without considering MHD turbulence, 
other possibilities remain open (c/GW06). Here we assume weak, but widespread turbulence 
throughout the asteroid formation region. 

Turbulence is an essentially lossless cascade of energy from large, slowly rotating ed- 
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Symbol 



Definition 



Equation or section 



l,v{l),te{l),oj{l) 


eddy scale, velocity, lifetime, and frequency 


sec. 2.2 


L,VL,nL 


largest eddy scale, velocity, and frequency 


sec. 2.2 




Kolmogorov (smallest) scale and lifetime 


sec. 2.2 


a, fl, Vk = 0,^ 


distance from sun, orbital frequency, Kepler velocity 


sec. 2.2 


H, c, pg 


gas vertical scale height, sound speed, and density 


sec. 2.2 


Re 


Reynolds number 


sec. 2.2 


a, ut 


kinematic viscosity Ut = acH 


sec. 2.2 


ts 


particle; stopping time 


eqn. 1 


Pp 


local mass density in particles 


sec. 2.3 


$ 


local mass loading factor — Pp/pg 
normalized gas enstrophy uj'^{l)/ 


sec. 2.3 


S 


sec. 2.3 


N 


cascade level corresponding to lengthscale I 


eqn. 2 


m,p{m) 


cascade multiplier and its PDF 


sec. 2.3 


P{^,S) 


joint PDF of mass loading and enstrophy 


sec. 2.3 




P($, S) at the peak of an IMF 


sec. 3.3 


Pgoal 


value of P* needed to create Mpa 


eqns. 10-15 




conversion timescale of mass into planetesimals 


sec. 3.3.1 




values of $,iV at P* 


sec. 3.3 


M 


mass accretion rate of gas 


sec. 3.3.1 


Mpa 


mass accretion rate of planetesimals 


eqn. 9 


A, Aq 


actual and canonical solids abundance relative to gas 


sec. 3.3 


tG 


dynamical collapse time of a dense clump 


eqn. 3 


^sed 


sedimentation time of a dense clump 


eqn. 4 


WeG,We*G 


Gravitational Weber number and its critical value 


sec. 3.1 


P 


pressure gradient parameter 


sec. 3.2 


a{a),pg{a),H{a),P{a) 


radially dependent nebula properties 


eqns. 5 




reference distance from sun (2.5 AU) 


sec. 3.2 


PR 


Roche density 


sec. 3.2.1 




thresholds for primary accretion 


eqns. 6-8 


Fv,Fp,Ft{> T) 


volume, particle, and time fractions exceeding threshold T 


sec. 3.5.1 


tenc 


particle encounter time with planetesimal-forming clump 


eqn. 17 


Aa 


radial diffusion (mixing) extent in tenc 


sec. 3.5.1 



Table 1: Symbols, parameters, and functions used in this paper 
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dies with lengthscale L and velocity Vl, which are forced by (currently unknown) nebula- 
scale processes, through smaller and smaller scales of size Z, having correspondingly shorter 
eddy timescales te(0' some minimum lengthscale 77, called the Kolmogorov scale, where 
molecular viscosity Vm can dissipate the macroscopic gas motions and turbulence ceases. 
We characterize the intensity of turbulence by the parameter a which sets the disk tur- 
bulent viscosity vt = LVl = acif, where c is the gas sound speed, H is the nebula ver- 
tical scale height, L — Ha^/'^, and Vl — ca}-!'^ . Then the turbulent Reynolds number 
Re — {L/r})^/^ — acH/v^. A typical T Tauri-hke nebula with mass accretion rate M ~ 
a fewxlO~^M0/yr, channeling 2-3% of its accretional energy into turbulence, would have 
a ~ 10-^ - 10-3 or Re = 10^ - 10^ at 3 AU (Cuzzi et al 2001; henceforth CHPDOl; also 
CW06). One may distinguish between turbulent viscosity and turbulent diffusivity (Prinn 
1990): the former is problematic in, for instance, convective turbulence (Ryu and Goodman 
1992) but the latter is robust in turbulence of all kinds, and it is the latter that drives our 
primary accretion scenario. A significant nebula turbulent diffusivity can also help explain 
the persistence of ancient, refractory inclusions in chondrites (Cuzzi et al 2003, 2005) and 
the abundance of crystalline, moderate volatility silicates in the STARDUST sample (Ciesla 
2009; cf. also Bokelee-Morvan et al 2002). That is, nebula turbulence can mix material 
radially by significant distances over time. 

In most cases of realistic, high- it!e turbulence, the Kolmogorov energy spectrum is a good 
approximation, where for a wide range of lengthscales r] < I < L, the turbulent kinetic energy 
density E{1) is given by the inertial range expression E{1) = {Vl/2L){l/L)~^/^. The eddy 
frequencies then scale as uj(l) = l/te{l) = v{l)/l = {2lE{l) f''^/l = VLl{1/ L)-^'^ , where v{l) 
is the velocity of an eddy of size Z, and the large eddy frequency is usually identified as 
the local orbit frequency (CHPDOl, Johansen et al 2007). These properties tend to be 
independent of the forcing mechanism and even of the Reynolds number of the turbTilcncc. 
Even if the initial forcing is anisotropic (as perhaps for MRI turbulence), smaller eddies 
become more isotropic as the 3D nonlinear cascade proceeds (Kato and Yoshizawa 1997). 
High Re, inertial range turbulence is sufficiently scale-free (Falkovich and Sreenivasan 2006) 
that using statistical and spectral properties from limited inertial ranges to characterize more 
extensive ones (those at higher Re) is an appealing approach. We make extensive practical 
use of this "cascade" property, as described in section 2.3. We note here for future use in 
section 3.3.1 that the Kolmogorov eddy timescale = l/{ilL{L/T])'^^^) = l/{ilRe^^'^). 

Particle-gas interactions: Particle interactions with the gas are characterized by the 
particle stopping time ts which, for particles of interest here, is defined by the Epstein drag 
law: 

ts = rps/cpg, (1) 

where r and ps are particle radius and internal density, and c and pg are the gas sound speed 
and density (see CW06 for more discussion). 

Particles interact with the gas, turbulent or not, within their stopping time and acquire 
inertial space (absolute) and random (relative) velocities accordingly. The relative velocities 
between particles determine the outcome of their collisions (sticking, erosion, or breakup), 
and the inertial space velocities determine the degree to which they diffuse radially and 
vertically, thus controlling their settling to the midplane (see eg DubruUe et al 1995, Wei- 
denschilling 1997, Ormel et al 2008, Brauer et al 2008). In the dense midplane layers of 
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cm-m size particles which can form in nonturbulent nebulae, the local gas is driven to near- 
Keplerian speeds and relative velocities between particles remain low enough for continued 
growth to planetesimal size to occur very rapidly, with or without the help of various mid- 
planc instabilities, on timcscalcs of 10'^ — 10^ years (Cuzzi et al 1993, Weidenschilling 1997, 
2000; Youdin and Goodman 2005). This is actually problematic in view of the extended 
formation epoch of primitive bodies discussed above - the process may go to completion too 
rapidly (section 2.1; CW06). 

Even weak turbulence, however, frustrates growth at some limiting size which depends on 
the local gas density (Dominik et al 2007). As particles grow they become more vulnerable 
to mutual destruction, because their increasing stopping time couples them to eddies of 
increasing size and velocity. In turbulence this coupling is captured by the Stokes number 
St — tgOU where lu can represent the eddy frequency on any scale - commonly either the 
large eddy scale L (Q^ ~ Q) or the Kolmogorov scale rj. Particles with tg comparable to 
the lifetime of the largest eddies (t^Q ~ 1) achieve the highest velocities Vl — a^^^c. For 
a ~ 10^^, such particles collide at relative velocities Vrei ~ V/, ~ 10'^ cm/scc - which are 
probably disruptive (Stewart and Leinhardt 2009). For a range of nebula properties, particles 
in the dm-m radius range have this property (see CW06, figure 1, or Ormel et al 2008) and 
this problem is commonly referred to as part of the "meter-size barrier" . However Brauer et 
al (2008), who assume a relatively low gas density, see growth frustrated at an even smaller 
size because the lower gas density leads to longer tg. Recent lab work is challenging some of 
the sticking assumptions of prior years at even lower velocities (Giittler et al 2010); models 
using these new results even suggest that a "bouncing barrier" might preclude growth beyond 
objects having masses not much larger than those of chondrule precursors (Zsom et al 2010) 
- again, depending on nebula properties. Obtaining and retaining objects with igfl ~ 1 is 
perhaps the major issue in the primary accretion scenario of Johansen et al (2007), which 
relics on an abundance of such particles in moderate-intensity turbulence because they drift 
rapidly into high-pressure ridges to become concentrated. Progress in this area will be 
interesting to follow. Meanwhile, we focus on a different accretion pathway, that also relies 
on turbulence but acts on particles much smaller than a meter, which are excited to small 
relative (collision) velocities well below the disruption threshold (Hogan and Cuzzi 2003, 
Ormel and Cuzzi 2007), and have sizes directly relevant to meteorites. 

2.3 Turbulent concentration and the cascade model 

Small particles diffuse in turbulence, but the trajectories of particles of a certain well-defined 
aerodynamic stopping time avoid fluid zones of high vorticity and converge in zones of low 
vorticity. Here, concentration factors C = Pp/p^ may be ^ 1, where Pp and p^ are the local 
and nebula-averaged particle mass density, respectively. We define the local mass loading 
$ = Pp/ Pg, where $ can thus also be ^ 1. The maximally concentrated particles have a 
stopping time tg equal to the overturn time t^^ of the smallest eddies (which have size 77, the 
Kolmogorov scale). Two "fingerprints" of this turbulent concentration (TC) seem evident 
in the meteorite record. The typical chondrule size (crudely, mm-diameter) is naturally 
explained by TC merely by requiring tg — tn (CHPDOl). An equally compelling fingerprint 
is the very characteristic chondrule size distribution, which is very similar across meteorite 
groups when scaled to the mean size, and is an excellent fit to the distribution predicted by 
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TC (CHPDOl). It has recently been shown that the TC size distribution is even statistically 
preferable to a lognormal distribution (Teitler et al 2009). It should be noted that, in the 
outer nebula where gas densities are lower and turbulent intensities plausibly larger, much 
smaller solid grains, or, more likely, porous aggregates of grains, such as seen in cometary 
IDPs, become the preferred candidates for TC rather than chondrules, for which tg would 
be too large (CHPDOl; see their section 3 and figure 1). Of course, it would be nearly 
impossible to extract "fingerprints" of the process in, for instance, returned KBO samples, 
after porous aggregates had become compacted in a parent body. 

Cascade model: The spatial distribution of $ = Pp/pg is determined only statistically, 
and must be studied with Probability Distribution Functions (PDFs) which depend on the 
nebula Reynolds number and the spatial scales of interest (CHPDOl). Because the nebula Re 
is far higher than achievable with current 3D fluid models, we have developed and employed 
a "cascade model" which, while not reproducing the physical structure of turbulence (vortex 
tubes and the like), has been shown to reproduce the PDFs of a number of attributes of 
turbulence (Mcnevaux and Sreenivasan 1991, Sreenivasan and Stolovitsky 1995). This model 
was described in detail by Hogan and Cuzzi (2007), and summarized by CHS08, so will be 
even more briefly sketched here. 

In turbulence, a number of properties (energy, velocity, vorticity, and particle abundance) 
can be thought of as being partitioned unequally and losslessly into sub-elements of eddies 
as they bifurcate. The partitioning fractions at each bifurcation are taken as m and 1 — m, 
where the "multipliers" m are drawn from a PDF p{rn) which is generally independent of 
eddy scale throughout the turbulent inertial range (Meneveau and Sreenivasan 1991, Juneja 
et al 1994, Sreenivasan and Stolovitsky 1995; see however Bee et al 2007 where some evidence 
is presented for scale-dependence in the context of preferential concentration). We determine 
the PDFs of these multipliers p{m) from our highest Re 3D models, which still cover only a 
limited range of eddy scales or bifurcation levels (Hogan and Cuzzi 2007). Each bifurcation 
is thought of as a level in a cascade; in the cascade model, we extend the multiplier process 
to even deeper levels (which one can think of as the smaller eddy scales achieved at higher 
Re). Unless ni ~ 0.5, repeated application of asymmetrical partition fractions (m, 1 — m) 
constantly creates more extreme values (higher and lower) of all parameters as the number 
of levels increases; this is referred to as intermittency - the local value becomes not more 
well-defined, but more highly variable at smaller scales (see the readable introduction by 
Meneveau and Sreenivasan 1991). For a cube, three orthogonal ID bifurcations, or levels, 
are needed to generate 8 subvolumes of linear size Ij^i = lj/2, and thus Re^^^ = L/rj — 
2^/3 = 10^°^^'^/^ ~ 10^/^*^, where is the total number of levels in the cascade. The general 
cascade relation giving the lengthscale associated with a given cascade level N , applied to a 
nebula situation with a large eddy scale L, is thus 

/ = 2-^/3^ = 2-^/3^^1/1 (2) 

Cascade models can achieve much higher Re than Direct Numerical Simulations (DNS); to 
match our full 3D DNS simulations at Re = 2000, the cascade model needs only about 15 
levels, taking about 10 cpu-hours (for 1024 realizations) compared to over 90000 cpu hours 
to converge a single full 3D simulation. 

Our particle-gas cascade model (Hogan and Cuzzi 2007) simultaneously treats $ and 
local enstrophy S — uj'^{l), where uj{l) is a vorticity on lengthscale I, using two distinct sets 
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Figure 2: The PDFs -P($, S) for four different values of Re, computed from 3D direct numerical simulations 
(solid contours) are compared with cascade model predictions for the appropriate number of levels (dashed 
lines; N=9, 12, 15, and 18 respectively, corresponding to the various values of Re). Note that, as Re and 
N increase (from panel (a) to (d)) the variance of both $ and S increases (more extreme values of both 
are seen). Mass loading has begun to affect the PDFs at the higher $ values and deeper cascades plotted. 
Figure from Hogan and Cuzzi (2007). 



of multipliers, and allows for the observed spatial anticorrelation of $ and S* on a statistical 
basis. The results of these cascades are binned into a second kind of (2D) PDF P{^,S). 
Examples are shown in figure 2. The meaning of P($, 5) is volume fraction (volume per 
unit nebula volume) having a particular combination of particle mass loading factor $ and 
relative enstrophy S — uj^{l)/ (w^(/)), where {uj'^{l)) is the average enstrophy at scale size 
/. P($, S) is given per unit logio($), P^r unit logio(5') and differs slightly in meaning from 
expressions in Hogan and Cuzzi (2007; see Appendix). P($, 5) is a function of level in 
the cascade, because going to deeper levels (smaller scales I) always enhances the variance 
of its properties (Meneveau and Sreenivasan 1991; figure 2). The two-dimensional nature 
of P($, S) is also essential; we will show that the threshold conditions allowing planetesimal 
formation depend on both $ and S, as well as level (section 3.2). 

As the particle mass density increases relative to the gas mass density, it affects the 
physics of turbulent concentration. Hogan and Cuzzi (2007) showed how particle mass load- 
ing affects the cascade; multipliers for mass loading $ and enstrophy S are shown to depend 
on the local mass loading itself. As mass loading increases towards $ ~ 100, multiplier 
PDFs p(m) narrow towards a delta-function at m = 0.5, implying an equal probability of 
partitioning and no further trend to intermittency (see the discussion in CHS08 or Hogan 
and Cuzzi 2007). In this situation, there can be no further increase of $ as the cascade 
level increases, and $ saturates near 100. The reason for this is not completely understood, 
but it is probably due to a combination of inertial effects (conservation of kinetic energy) 
and dissipation of turbulent kinetic energy (which may be smaller than in terrestrial ex- 
periments because of the tighter coupling of the particles to the gas in our regime). The 
mass-loaded cascade models of Hogan and Cuzzi (2007) showed good agreement with actual 
3D two-phase, mass-loaded, DNS results (figure 2). 
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An important property of the cascade models is that the cascade level N can be related 
to the corresponding nebula lengthscale I, for some given nebula a (equation 2 above), and 
we can thus calculate the mass of the planetesimal which forms from a given volume of 
material at each level. In section 3 we will show how the cascade model, combined with 
certain thresholds for the stability of dense clumps, leads to the primary mass distribution 
function, or IMF, of planetesimals. The results noted in section 3.1 (see also CHS08) point 
to spatial scales of interest which are 10^ — 10^ times larger than the Kolmogorov scale rj 
emphasized by CHPDOl; that is, because of the role of mass loading, we no longer beheve 
the theoretically high concentrations of CHPDOl are achievable at small scales, and thus, for 
reasons discussed in section 3.2.3, Icngthscales as small as 1] also become less relevant. There 
are fewer cascade samples in shallower cascades - larger / means lower N and the number of 
samples is about 2^ for a cascade with N levels or bifurcations. Thus, we needed to run many 
cases to obtain the proper statistics to explore the low-P ranges of ($, S) that exceeded our 
newly defined thresholds (section 3.2) and were capable of becoming planetesimals. We ran 
cascade models for a period of several months on NASA's HEC Altix and Origins computers 
at Ames, ultimately running 10^, 10^, and 10'' cases at 24, 20, and 15 levels to build up 
statistics such as shown in figure 3 (next section). 

3 Determination of Initial Mass Functions 

We first review how the self-gravity of a clump enters, which is more subtle than usually 
believed. We then derive three different thresholds on different combinations of clump density 
and size, and local vorticity, that determine which dense clumps can become primary sandpile 
planetesimals in a turbulent nebula environment. 

3.1 The role of self gravity 

Gravitational Instability (Gl) or inexorable collapse on a dynamical timescale 



is a well-used tool in the cosmogonist's toolbox, but we have found that the traditional 
concept of Gl is not appropriate for small particles which have stopping times tg much 
less than dynamical times to- CHS08 recently rediscovered numerically a result originally 
obtained analytically by Sekiya (1983), and since apparently forgotten: that gas pressure 
stabilizes dense clumps of particles against traditional gravitational instability on dynamical 
timescales (Gl). In the regime where particle-gas coupling is strong, particles which begin to 
collapse under their self-gravity drag and compress the entrained gas, producing a radial gas 
density and pressure gradient, which in turn prevents the gas and tightly coupled particles 
from undergoing Gl until the particle mass loading is 10^ times larger than the traditional Gl 
criterion. Sekiya (1983) called the mode of particles and entrained gas that arises under these 
conditions, ordinarily assumed for traditional Gl, a 3D "incompressible mode" of instability. 
Within such blobs, Sekiya suggested and CHS08 showed that particles of radius r can only 
sediment slowly inwards at their terminal velocities, on the timescale 



to = (4G$p,) 



1/2 



(3) 



tsed = l/4:G^pgts = C 



■./AG^rp, 



(4) 
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on the order of 10^ — 10^ orbit periods for typical chondrules and $ = 100 — 10. 

CHS08 explored the ability of such dense clumps to resist disruptive forces for a time 
on the order of tsed- As a clump settles vertically towards the midplane under the vertical 
component of solar gravity, or orbits at a velocity near Kcplerian {Vk), it incurs a ram 
pressure from the nebula gas. In the latter case the headwind arises because the gas, being 
pressure-supported, orbits more slowly than Keplerian at a speed of (1 — (3)Vk where (3 ~ 
10~^ (see next section; Weidenschilling 1977, Nakagawa et al 1986; CHS08). Conservatively 
assuming each clump moves as a Keplerian object, the mean settling speed for a clump 
formed at altitude z above the midplane is Vz ~ {z/ajVx ~ {z/ H)f3^/'^VK since (3 ~ {H/aY- 
Then requiring the vertical headwind be smaller than the azimuthal headwind implies z/ H < 
Within a vertical distance Hf3^^'^ of the midplane, the vertical settling velocities are 
small compared to the orbital velocity difference between the pressure-supported gas, and 
the azimuthal ram pressure on a strengthless Keplerian clump dominates all other disruptive 
forces (CHS08). In section 3.3.1, we will restrict the volume in which plantesimals can form 
to this near-midplane region. 

CHS08 developed a toy model based on an analogy with the Weber number We in the 
familiar raindrop problem, where We is the ratio of surface tension to ram pressure forces 
acting on a fluid droplet moving at velocity AV relative to a less dense fluid. They defined 
a "gravitational Wcbcr number" Wea which balances the ram pressure force per unit area 
with the self gravitational force per unit area of a strengthless clump of initial diameter / and 
particle density Pp. The premise was that certain combinations of diameter / and particle 
mass density Pp = ^Pg would stabilize a clump against being disrupted by a headwind of 
magnitude AV — ^Vk- CHS08 ran a range of numerical models of clumps experiencing a 
steady nebula headwind from the more slowly orbiting gas, to validate the toy model, and 
determined that stability was indeed achieved for Wee greater than some critical value We*Q 
of order unity. Viscous losses of material around the periphery of their numerical clumps 
limited their numerical runs, but such large viscous erosion is an artifact of the numerics 
and would not be present in the actual nebula case. They noted that the combination of 
parameters required for stability of a dense clump implied a substantial size for the ensuing 
sandpile (10- 100km radius), and pointed out the similarity of this size to the "fossil asteroid 
belt" modal size of Bottke et al (2005). Below we show simplistically, but quantitatively, 
how a combination of dense clump stability thresholds may determine the IMF of primary 
planetesimals. 

3.2 Thresholds for primary accretion in (5, $) space 

Our prediction of primary object IMFs is based on mapping three different kinds of threshold 
onto the cascade probability contours (figures 2 and 3) . These are not thresholds at which 
any sort of traditional "fast" instability occurs (section 3.1) - rather they are thresholds 
which allow dense particle clumps to avoid ram pressure disruption by the nebula gas for the 
long time {tged ~ 10^ — 10^ orbits) required for the particles in them to sediment into their 
mutual center, creating a sandpile planetesimal. The key step in deriving planetesimal IMFs 
is connecting the thresholds derived below (functions of lengthscale I) to the cascade model 
PDFs (functions of level A^), which we do using the cascade relation given in equation 2. 
In the following sections we incorporate simple scaling of our criteria with distance a from 
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the sun, based on powerlaw approximations for the nebula gas surface mass density a{a) = 
(To{a/ao)~^ and, to a less important extent, mean temperature T(a) = To{a/ ao)~^ ■ The 
combination of radial density and temperature gradients leads to a generally outward radial 
pressure gradient dP/da which, normalized by the gravitational coriolis force, is specified 
by the nondimensional parameter /3 = {dP / da) / {2pgaQ?) (previous section, Weidenschilling 
1977, Nakagawa et al 1986; possible complications are discussed in section 3.4). Then, 
following Cuzzi et al (1993, equations 54-59): 

a{a) = 2H{a)pg{a) (5) 
Pg{a) = p,(ao)(aK)-(2p-'/+3)/2 ^ p^(a,)(a/a„)-(^'+3/2) 

H{a) = i/(a<,)(a/a„)(3-'?)/2 ^ H{ao){a/aof^ 

fi{a)^fi{ao){a/aof-'i 

where the g-dependence is weak for g ~ 1/2 and ignored for simplicity except in ^; it can 
be easily allowed in more detailed studies. We also adopt Vt[a) = il{ao){a/ ao)~^^'^ , and let 
flo, Po, Ho, Pgo — i^(ao), /3(oo), etc. take their nominal values at ao—2.5A\J. 

3.2.1 Threshold $i: Rotation and gravitational binding 

The first question most people have is, are the clumps rotating too quickly to be bound? 
This threshold is determined by comparing the local gravitational timescale to (equation 
3) and the local eddy timescale l/u!{l), where the local eddy frequency u{l) is treated as 
a vorticity. This threshold dominates when the local vorticity exceeds the global value 
(see section 3.2.2). Requiring tc < 1/^(1) (Toomre 1964, Goldreich and Ward 1973) is 
conservative here because eddies don't truly "rotate" (many times) with timescale l/uj{l); 
rather, l/uj{l) is their existence lifetime before bifurcating. Even though the dense zones 
of small particles of interest here cannot collapse on the timescale to (CHS08; section 3.1), 
they can become bound entities based on a criterion close to this (Sekiya 1983). Then 
to = (4G$pp)"^/2 < l/id{l) or $ > uj'^{l)/AGpg determines our first threshold $i. To express 
$1 in terms of S* = oo^{l)/ (a;^(/)) we use the inertial range mean enstrophy on scale /, 
(a;^(/)) = Vl\[L/IY/'^ (section 2.2) where VLl is the large eddy frequency, generally taken to 
be the orbit frequency Vt. We then use = GMq/o^ where Mq is the Sun's mass and a 
is the distance from the Sun, and also the definition of the Roche density pr = 3M0/47ra^ 
(Safronov 1991) to get $i(5) = {pn/ Pg){L/ 1)^/^3. This relat ion is extended to arbitrary 
semimajor axes a using the a-dependcncc of pR^a) / pg(a). Then Pr/ Pg = Kq^q,/ aoY~^^'^ , 
where Kq = (ZMQ/AnpgoO?^). We then use the cascade relation / = 2~^/^L to express 
{L/l)^/'^ — 2^^/^. In the nebula, L — B.a^l'^ is the large eddy scale. Combining these 
relations leads to 

^r{S,a)^2'''l'K,s[^^ . (6) 

Note above that $i has no explicit a-dependence, but each iV implies a lengthscale / which 
does depend on a. Each threshold $i appears as a diagonal line in figure 3, colored 
according to its value of N . 
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The more refined stability analysis of Sekiya (1983) is easily generalized to this situation. 
Sekiya assumes the relevant rotation frequency is the orbital frequency Q and finds that 3D 
incompressible modes become marginally bound at AnGpp/Q'^ ~ 10 (his section 4 and figure 
3). We generalize this to AnCpp/ = 47rG^Pg/u'^{l) > 10, and follow the same logic as 
above, substituting = {u'^{l)) S = VL^{L/lY/^S = and relating il^ to pn/pg 

as above, to obtain a threshold value ^sek = t^I' Because even this more sophisticated 
analysis is idealized itself to some degree, we explore the imphcations of this factor of 10/3 
as one example of the uncertainty of the predictions in the figures and tables of section 3.4 
below. 



3.2.2 Threshold Smin- global rotation 

Notice in figure 3 that the cascade PDFs extend to very low values of relative enstrophy S ~ 
fewxlO"^. However, on the long timescales tged a clump cannot be guaranteed of remaining 
in fluid zones with such low vorticity, and will experience the global rotation as a minimum. 
We thus impose a minimum local vorticity given by the global shear rate Q. We express this 
in terms of S as 

S = > Srmn = = n'/{2'^"^,f = 2"^^/^ (7) 

Thresholds of Smin appear as vertical lines in figure 3, colored according to N . Regions lying 
to their left can be disregarded as candidates, having unrealistically low S to characterize 
the long timescales involved in sedimentation to sandpiles. 



3.2.3 Threshold $2: Ram pressure and the Gravitational Weber Number 

As described in section 3.1 and CHS08, self-gravity of a dense clump can play the role 
of surface tension and stabilize a clump against the disruptive ram pressure of the nebula 
headwind if $2^ > l^aVt/ {2GpgWe}.Y''^ (CHS08 equation 4). CHS08 suggest that the most 
favorable region for clump survival is within some small vertical distance jS^^^H of the mid- 
plane, where vertical settling of dense clumps under solar gravity is negligible and only the 
azimuthal headwind remains (section 3.1). The limited, coarsely gridded numerical simula- 
tions in CHS08 were unable to establish a precise value for Wcq, but it appears to be of 
order unity, which we adopt here; for raindrops falling in Earth's atmosphere, Wcq — 8, so 
refining this constant in the nebula application is worthy of more attention. As above we sub- 
stitute fi^ = GMq/o,^, and, closely approximating the definition of pr as M0/(4a^), obtain 
$2(0)/ = /3a{2pR/pgWe*Gy/^. Then using pn/pg from section 3.2.1, I = from the 

cascade relation, and scaling a/H and ^ with a as above, we obtain in a straightforward way 



^ Ho J \aWeQ J \ac ^ 

where Kq is defined in section 3.2.1 and again we retain the A^-dependence. Note that $2 is 
not a function of S, and thus appears as a horizontal line in figure 3 for each value of A^, but 
is an exphcit function of a because of the I factor in the threshold equation for $2^ (CHS08 
equation 4, and above). CHS08 discuss why other possible gas effects, such as turbulent 
pressure fluctuations, are negligible compared to simple ram pressure. 
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3.3 Derivation of the Initial Mass Functions 

The two most important things about an IMF are (a) the shape of the mass distribution 
P{M), in particular its modal value if any, and (b) its absolute value, giving the rate at 
which primary planetesimals of those masses are created. Our cascade model is the key to 
both. The cascade model PDFs P($,S') (section 2.3) refer to scale / = where for 

the nebula L ^ Ha^^"^ is the largest eddy diameter. Thus the cascade model provides us 
with the volume density (the occurrence probability or volume per unit volume) of zones 
having a particular combination of density and vorticity on a specific nebula lengthscale /, 
and allows us to calculate both a mass associated with each bin (given by M = $pg/^), and 
the abundance of these bins at any time (given by P($, S*)). As we increase the cascade level 
A^, we sample statistics at smaller /, where the PDF is more intermittent and the probability 
contours expand, limited by the constraint of saturation near $ ~ 100 (section 2.3; Hogan 
and Cuzzi 2007; figure 3). The values of P($, S) depend explicitly on N and the initial value 
of the total solid/gas ratio A, with canonical value Ao which we take to be 10~^ everywhere 
for reference, comprising particles with sizes suitable for turbulent concentration [tg ~ t,,; 
section 2.2). 

The thresholds $i, $2, and Smin (equations 6-8) also increase with N (a result of their 
/-dependence). For any N, the most common planetesimal mass is that corresponding to the 
clump (of size I) with the highest value of P($, S) = lying along the threshold lines $1, 
$2, for S > Smin- As A^ changes, the contours and thresholds evolve at different rates; thus 
P/v varies with A^ and there is typically some maximum P/v = P* at some value of A^ = A^*. 
This defines the peak of the distribution at A^*, $*, P*. Because there is also a mass M 
associated with any ^,1{N), and pg, Pn{M) provides the complete IMF and has a modal 
mass M($*). 

This situation is best perceived in a sequence of snapshots at different A^, which are 
difficult to present in the format of a printed page (see online supporting material or 
|http : // spacescience . arc. nasa .gov / media /st aff /j eff-cuzzi / IMF .ppt ) . We attempt to present 
it in figure 3 using two different colors for the contours and thresholds associated with two 
different values of A^. The planetesimal diameters plotted in the right panel are derived from 
M = ^*pgl^, assuming a planetesimal density of 2.0 g/cm^. 

To summarize, the placement of the thresholds $1, $2, and Smin depends on A^ and the 
physical parameters of the nebula model assumed: the nebula a, the local gas density pg or 
surface density a, and the headwind parameter p. The placement of the contours P($, S) 
depends on A^ and the local solid/gas ratio A, which can be enhanced over cosmic abundance 
Ao (here assumed to be 0.01). Primary IMFs vary accordingly. Extension of the theory to 
30 AU is straightforward (section 3.2 and 3.3.1), depending on the radial dependence of a{a) 
and /3(a), as determined by powerlaw relationships plus whatever (ill-constrained) radial 
variation there might be of turbulent intensity a. In section 3.4 we show preliminary IMFs 
at 2.5 and 30 AU, for a number of nebula parameters. Before describing these, we outline 
our approach to constraining the vertical component of the IMFs - the actual creation rates 
as functions of size. 
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Figure 3: A closer look at how the IMF is determmed from the PDFs S) and the thresholds for 
primary accretion (sections 3.2 and 3.3; see also figure 2 and section 2.3). Left: Contours show the cascade 
model results for (the log of) fractional volume P{^, S), running from 10^^ to 10"^*^, at two different cascade 
levels N; the blue contours are for the larger N. This cascade refers to a case with A = lOAo- Also shown 
for the same two levels are the thresholds $i = ^sek (diagonal), $2 (horizontal), and Smin (vertical; see 
section 3.2). Note that, as the contours expand with increasing N, the thresholds recede up and to the left. 
Right: the curve in the right panel (section 3.3) plots the maximum value of -P($, 5") in the region of ($, S) 
that exceeds all three thresholds, as a function of calculated diameter D, with one point for each value of 
N; red and blue points refer to the same cases in the left panel. The curve traced out as N varies gives the 
primary accretion IMF, which has a peak at some N = N* , defining P* and the associated The curve 
in the right panel can be associated with the case pgo = 3 x 10~^, a = 10"'*, /3o = 10"'' (figure 4b); as noted 
in section 3.4, these values differ from canonical values. 



3.3.1 Planetesimal creation rate 

A successful model must reproduce estimates of the mass originally created in primordial 
planetesimals in some region, over the time available; this is the average primary accretion 
rate Mpa- We will compare Mpa from our models with expectations for the solar system. 
Primary accretion occurs when clumps having local density of solids $Pg, occupying some 
small volume fraction of the nebula P{^,S), become stable against disruption and form 
a cohesive planetesimal in their sedimentation timescale tsed = ^/ ^G^pgts (section 3.1; 
CHS08). We can assume primary accretion is dominated by a region near the peak of 
each modeled distribution (a more refined approach is described in the Appendix). Then we 
require the set of parameters (P*, N*) at each peak or modal value to satisfy a stipulated 
primary accretion rate Mpa and solve for the value of P* which we refer to as the "goal" value 
Pgoai- We can then normalize the various distribution peak values P* by the corresponding 
Pgoai for the same parameter set, to assess how well the parameter set achieves the stipulated 
Mpa- The available nebula volume between semimajor axes ai and 02 is vr(a2 — a\) ■ 2i7/3^/^, 
where only some vertical fraction may be suitable for this process (section 3.1). Thus 

Mpa = ($*Pg)P*(27r(a2 - al)H 13^'^) /tpa- (9) 
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In equation (9), the numerator represents the total amount of mass at any instant lying 
in zones which exceed the various thresholds $2, 'S'^in) and can become planetesimals. 
The demominator tpa is the timescale on which primary accretion converts this mass into 
planetesimals. Wc then solve equation (9) for the values of P* = Pgoai which are needed to 
produce the estimated primary accretion rate as 



-J ^patpa ^pa^sed 



(10) 



where we thereby define a primary accretion rate in terms of surface density: dpa- 

In equation (10), for specificity and to be conservative, we adopt numerical values of 
tpa = tsed for the relevant mass production timescale; this is not a well-defined selection but 
assumes that all the physics of clump formation and dispersal, including various dynamical 
and fiuid timescales, is captured by the ensuing average volume fractions P($, S), such that 
the rate at which sandpile planetesimals appear is then simply the proto-sandpile mass so 
defined at any given time, divided by the time it takes them to become sandpiles. Because 
this choice is uncertain by a large factor, we will carry a final factor of tpa/tged which will 
illustrate the sensitivity of our results to the uncertainty in tpa- While studying the results 
presented in the next section, the reader should keep in mind that the value of Pgoai would 
be significantly smaller, and thus the normalized IMFs in figures 4-6 would be considerably 
closer to unity, if the timescale in the denominator of equation (9) were, instead of tgedi 
the formation time of a clump (plausibly on the order oft/, ~ I/^^l or roughly the orbit 
time, which is a factor of 10^ — 10^ shorter than tged)- Indeed Chambers (2010) has assumed 
an even smaller timescale for tpa-, comparable to the (shorter) eddy timescale at lengthscale 
I « L. The question of the most appropriate approach to estimating Mpa is a fruitful 
subject for future consideration. 

Below, we will use crude estimates of Opa at 2.5 and 30 AU to constrain our model 
predictions of P{^,S). First, we rewrite equation (10) in a more useable form, combining 
all occurrences of familiar nebula parameters. Wc assume the particle stopping time for 
preferential concentration t^, which occurs in the definition of t^ed, is equal to the Kolmogorov 
eddy timescale tr, which depends on nebula properties (see section 2.3): 

t ^t = ^ = = ^ (11) 

where Ki = 5.3 x 10^^ g-*^/^ cm~^, and we have expressed the gas kinematic viscosity as 
= ^H-z'^-/ (^H'lPg whcrc and ct/Zj are the mass and cross section of a hydrogen molecule, 
respectively (Cuzzi et al 1993). Chondrule-like particles satisfy this relationship in the 
asteroid belt region, but in the outer nebula, much smaller or less dense particles will be 
optimally selected (CHPDOl). Then substituting for tg we obtain 

&pan{apgHy/^ ^ &pana^/^ / tpg \ 

We now scale all radially variable quantities assuming nominal powerlaw nebula surface 
density and mean temperature discussed above, and obtain 

p _ ( \ f ^paa'/\a) \ fa\ ( t_pa\ 
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Regarding apa, a consensus belief is that the 2-4AU region of the primordial asteroid 
belt (prior to dynamical clearing) contained planetesimals with a mass of about 2Mq (Petit 
et al 2001, Chambers 2004, and personal communication 2009), which isotopic age dating 
suggests formed over about 2Myr. The Kuiper belt is less well constrained but required 
perhaps 40 in planetesimals between 16-30AU (Tsiganis et al 2005; sec sections 2.1 
and 3.4.2); for our scenario to be relevant this also must have happened before the gas 
vanished. Standard nebula hfetimes of about 3Myr (Haisch et al 2001) refer to the presence 
of warm dust, probably more relevant to the asteroid belt region than the Kuiper belt region. 
Currently, the lifetime of outer nebula dust can only be limited crudely to less than 10-30Myr 
(Carpenter et al 2005). For simplicity here, we simply assume the same accretionary lifetime 
for the outer nebula as for the inner nebula (2Myr); this is shorter than found by traditional 
incremental growth models (see section 2.1). Then (jpa(2.5AU) ~ 2.4 x 10~^^ g cm~^ sec~^, 
and (7pa(30AU) ~ 10~^^ g cm~^ sec~^. After some algebra, equation (13) becomes 

At 30 AU the results depend on the radial scaling parameters p and q (now embedded in K2 
below); we assume q — 0.5: 

..,30... = ..)(^)-(^)X3lB"^(£) 

where K2{p) = 10-5(a/ao)(3p-0-5)/2. thus 7^2(0.5) = 3.6 x 10-^ K2{1.0) = 2.4 x 10~^, and 
K2{1.5) = 1.5 X 10^^, and we have suggested (different) plausible values of a (sections 2.1 
and 2.2) at 2.5 and 30 AU, and a typical overall value of $*, for scaling purposes. 

3.4 Results 

Figure 4 shows preliminary IMFs we have derived at 2.5 AU, for a range of nebula pa- 
rameters, based on the methods described in section 3.3. The modal diameters for primary 
planetesimals fall within the 20-200 km range of uncertainty spanned by models of subsequent 
stages of evolution leading to the observed asteroids (Morbidclli et al 2009a, Weidenschilling 
2009; section 2.1 and figure 1). Each IMF is normalized by the "goal" value Pgoai for the 
combination of parameters defining each curve, calculated using equation (14) or (15). If 
the peak of the normalized IMF approaches unity, it implies that the case is capable of pro- 
ducing enough mass in planetesimals, in the time available, to satisfy current expectations. 
The actual values of P* and Pgoai ai'c tabulated in Tables 2-4. Clearly, some cases are more 
successful than others in this regard, but it is intriguing that the model even comes close 
to satisfying both of these independent constraints at once. However, to do so, the results 
shown for 2.5AU prefer a local background solid mass enhancement over cosmic abundance 
A/Ao — 10 (figures 4a, 4b), or A/Ao — 30 (figure 4c), and a headwind parameter /3 which is 
as much as 10 times lower than normally assumed {/3 ~ 10~^; Nakagawa et al 1986; Cuzzi 
et al 1993). Enhancement of solids over canonical values, and suppression of the headwind 
speed below canonical values, are not only important, they are connected (see below). 
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Comparison of figures 4a-c is instructive regarding the effects of uncertainty in the 
models. All results are obtained applying all three thresholds: $i, $2, and Smin- However, 
figure 4a shows results derived assuming the simple derivation of $1 in section 3.2.1, with 
A = lOAo- Figure 4b shows the implications of adopting, instead, the more refined threshold 
^sek = (10/3)$i (sec section 3.2.1) - the values of P* and Pgoai (and of the resulting nor- 
malized IMFs) decrease, while the modal sizes shrink slightly. In figure 4c, we show that 
increasing A/Ao by only a factor of three increases the IMF P* values dramatically, to the 
point that the normalized IMFs routinely exceed unity; thus for our assumption of tpa — tged 
(section 3.3.1), a degree of solids enhancement in the range A/ Ac ~ 10 — 30 is apparently 
called for. On the other hand, the normalized IMFs would increase by several orders of 
magnitude if we were to adopt tpa ^ instead of tpa ~ tged (see section 3.3.1 and Chambers 
2010). In section 3.5 we discuss the ability of these cases to match other constraints, where 
tpa plays no apparent role. 

3.4.1 Enhancement of local solids by radial and vertical decoupling 

Suggestions that A > Ao, where Ao is the cosmic abundance, are not new in the context 
of primary accretion. Growth to dm-or-m size (but perhaps no further) may be robust in 
turbulence (Dominik et al 2007, Ormel et al 2008), especially in the outer solar system 
where water ice might increase particle "stickiness" . This will cause large amounts of mass 
to migrate from the outer solar system to the inner solar system much faster than the gas 
evolves, elevating the relative abundance of solids significantly (Stepinski and Valageas 1996, 
1997, Cuzzi and Zahnlc 2004, Ciesla and Cuzzi 2006, Kornet et al 2001). A similar process 
was advocated by Youdin and Chiang (2004) in a nonturbulent nebula without particle 
growth. More recent models (Zsom et al 2010) which account for experimental results for 
silicates (Giittler et al 2010) find growth being frustrated at even smaller sizes, where radial 
drift would be much smaller - keeping material around longer in the inner solar system where 
ices are mostly absent. A complementary process is vertical settling of clumps formed at 
high altitude, which brings material to lower altitudes faster than would otherwise be the 
case (Wang and Maxey 1993, Aliseda et al 2002, Bosse et al 2006). We are actively studying 
both of these processes, and feel that a combination of them could lead to background 
abundance of A/Ag ~ 10 — 30 within z < where we have suggested primary accretion 

operates. Moreover, settling of solids towards the midplane can affect the local gas orbital 
velocity and thus the headwind experienced by particles (Nakagawa et al 1986). For plausible 
enhancements {A/Ag < 30) this is unlikely to lead to more than an order unity effect; solving 
for the headwind velocity using equations in Nakagawa et al (1986), for particles with short 
stopping times such as those of interest, gives a reduction in effective ;5 by a factor 1— p^/pg, 
or 0.7 for A/A^ = 30. Haghigipour and Boss (2003) showed in principle, and Johanscn ct al 
(2007) found in realistic 2D and 3D turbulent simulations, that the headwind parameter P 
can essentially vanish in local pressure gradient reversals, which may be long-lived but might 
only occupy a small fractional volume. Testing these speculations with actual models is an 
important goal for future work. 

The 30AU cases generally fall well below their respective goals if A — Ag is assumed 
(figure 5), especially for the steeper radial density distributions. However, some enhance- 
ment of solids at 30AU is not obviously out of the question, at least near the midplane: if 
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Figure 4: Initial Mass Functions (IMFs) at 2.5AU. The vertical axis plots the normalized function 
\og{P{^,S)/ Pgoai) (section 3.3.1), which peaks at the mode value P{^,S) = P* (section 2.3; right panel 
of figure 3). Peak values near unity indicate the scenario can produce the expected pre-depletion mass in 
planetesimals in the expected time. Figures 4a and 4b (top left and right) assume a background enhancement 
of sohds over cosmic abundance of A/Ao = 10, and figure 4c (bottom) assumes A/Ao — 30. Blue: a = 10~*, 
Red: a = 10~^, and Green: a = 10"'^. Figures 4b and 4c (and 5 and 6) increase our simply derived value 
of $1 (equation 6) by a factor of 10/3 to align it with the result of Sekiya (1983). The normalizing gas den- 
sities pgo refer to Oo =2.5AU. Curves are also labeled by the headwind parameter /3o(ao=2.5AU). Roughly 
speaking, larger a produces larger planetesimals, and increasing A/Aq or pgo and/or decreasing tpa or /3o, 
leads to higher production rates. These cases are tabulated in Tables 2-5. 
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Figure 5: Initial Mass Functions (IMFs) at 30AU, similar to figure 4; these assume different nebula 
radial surface gas density powerlaws a{ao){a/ao)~'' with p = 0.5, 1.0, and 1.5. The vertical axis plots 
\og{P{^, S) / Pgoai) (section 3.3.1). The curves here assume normal cosmic abundance {A = Ao) and the 
Sekiya value of $i. The normalizing gas densities pgo refer to 2.5AU, and curves are also labeled by the 
headwind parameter /3o(2.5AU). As in figure 4, larger a produces larger planetesimals, and increasing pgo 
and/or decreasing /3o leads to higher production rates. The scant results for p = 1.5 suggest we would 
need to extend the cascade model to higher N (smaller /) than we have so far, to capture the mode of the 
distribution, which would remain at low P*. Allowing for a smaller tpa could raise these normalized IMFs 
by a factor of 10^ — 10^ by decreasing Pgoai (section 3.3.1). 
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Figure 6: Same as figure 5; normalized IMFs at 30AU, assuming the Sekiya-adjusted value of $i and 
different surface density profiles -p =0.5, 1.0, and 1.5. These results differ from figure 5 by assuming a local 
abundance of solids enhanced over cosmic by a factor of 10 (A = 10 Aq). Normalized IMFs approach (or 
even exceed) unity for flatter radial distributions, and all normalized IMFs can be increased by decreasing 
tpa (section 3.3.1). 
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dense clumps are forming and settling, there will be some enhancement due to this process 
alone, as discussed above. Moreover, for nebulae that begin with >100 AU radial extent, 
some enhancement may occur due to particle growth at radii larger than 30AU, and sub- 
sequent inward drift, depending on the growth rate in these rarificd regions. Stepinski and 
Valageas (1997), Weidenschilling (2004), and Kornct ct al (2001) find enhancements of solids 
by factors of order unity by 0.1-lMyr, depending on a. On the other hand, Garaud (2007), 
using different assumptions (very low particle sticking), finds a depletion of solids at these 
distances. Further studies of the growth of solids in the remote outer nebula are important to 
continue because of these discrepancies and their implications for primary accretion. Over- 
all, we feel it is not implausible that some degree of enhancement of solids did indeed occur 
even at 30AU, although probably not as substantial as that which might have occurred at 
3AU. Because our results are so strongly dependent on the abundance of solids (figure 4c), 
we ran the same set of cases at 30AU assuming A — lOAg (figure 6). Indeed the normalized 
IMFs now come closer to, or even exceed, unity. It does appear that moderately steep radial 
nebula profiles, such as that of the traditional minimum mass nebula {p = 1.5) are seriously 
challenged to match the mass production rate goals as adopted here. Again, we note that 
smaller values of tpa also have the potential to increase the normalized IMFs closer to unity 
by several orders of magnitude (section 3.3.1; Chambers 2010). It is interesting to note that, 
in figure 5c, the normalized IMFs actually increase with increasing /3, so for the larger 
values of a ~ IQ-^ that may well apply to the outer solar system, nominal values of /3 might 
be acceptable (somewhat in contrast to the cases at 2.5 AU). 

3.4.2 Primary accretion efficiency and its impfications 

It is interesting to note that, at both 2.5 and 30 AU, only a small amount of the mass initially 
present actually gets accreted into planetesimals in this scenario; that is, our accretionary 
process is highly inefficient - quite distinct in this regard from the traditional 100% efficient 
minimum mass nebula assumption. Consider a typical p = 1 case with pgo = 10~^ g cm~^ at 
2.5 AU. Say our models achieving P* ~ Pgoai had A ~ lOAo on average. If this is a product 
of both radial and vertical settling, the solids surface mass density might be increased by 
a factor of 3 over our nominal nebula models (by, eg., inward radial drift of particles from 
further out). This would result in available mass of about I2OM0 in sohds in the 2-4 AU region 
and about QOOM^ in the 16-30 AU region. Then, P* ~ Pgoah or accretion of 2M0 in the 
2-4 AU region, means that only ~1.5% of the solids there were captured into planetesimals, 
and accretion of 4OM0 in the 16-30 AU region means that only ~4% of the solids there were 
captured. The balance (that is, almost everything) presumably escapes into the sun (or in 
the case of the outer nebula, escapes the solar system or helps feed the inner nebula) . 

There are some especially interesting implications for the Kuiper Belt. These efficiencies, 
while small, are still larger than those of traditional models of incremental accretion of the 
current Kuiper Belt {eg. Stern and Colwell 1997, Kenyon 2002). In these models, tens of 
M® of solids initially present in the region produced only the observed O.Ol-.lM® in KBOs 
- having an efficiency an order of magnitude smaller than that derived above because of the 
broad size distribution in which they are required to grow. Recent models of subsequent 
dynamical depletion (during evolution of Neptune) of the planetesimals that did form do not 
change this efficiency, but just require a larger starting mass. On the other hand, perhaps the 
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biggest remaining problem in the "dynamical emplacement" theory of the KBOs (Levison 
et al 2008) is the cold classical population, which has lower eccentricity, and perhaps higher 
binary fraction, than might be expected from dynamical emplacement from within 30AU. 
That is to say, growing at least the cold classical KBO population in situ, while emplacing 
the rest, might have some appeal; we suggest below that this is not out of the question given 
an edge in the solids density near 30 AU. These issues, and others related to a large initial 
mass in the Kuiper belt region, are described in section 2.1. 

As mentioned above, our primary accretion rates (and masses) are strongly dependent on 
A/Aq. Comparing figures 5 and 6 shows that a drop in A/ Aohy a. factor of 10 leads to a drop 
in P* / Pgoai by a factor of about 1000. Suppose there were a factor of 10 drop, or "edge", in 
the abundance of solids at around 30 AU, and all other parameters were slowly varying {eg., 
Weidenschilling 1997, Stepinski and Valageas 1996, 1997). We would then expect to form, 
over the 30-44AU region, 1000 times less mass in primary bodies than our goal value of AQM(^ 
- about O.O4M0. Indeed this is about four times the mass of the cold classical population 
(Bernstein et al 2004). Even if the dynamical evolution of Neptune removes about 90% of 
originally formed planctcsimals (A. Morbidelli, personal communication 2010), the scenario 
described in this paper is only a factor of 2-3 away from forming the required ~ O.IM® cold 
classical (primordial) population in situ. Moreover, this formation would not leave behind a 
massive population of other objects that would need to be eroded away and/or cause Neptune 
to migrate further than observed, because the bulk of the planetesimal mass formed by this 
scenario "forms big" at the mode in the mass distribution. Given the simplicity of the 
scenario presented here, and its sensitivity to uncertain parameters, it seems this possibility 
may be worth further thought. 

In future, more detailed studies, allowance should be made for the fact that a primary 
accretion mechanism designed to produce (part or all of) the current crop of KBOs in 16-30 
AU might also need to produce the roughly equal mass in planetesimals that grew into the 
cores of Uranus and Neptune, which would increase the mass accretion rate by a factor of 2 
or so; this does not strike us as prohibitive given the observational and model uncertainties 
involved. On the other hand, if Uranus and Neptune formed inside of 17AU {eg. as in 
Tsiganis et al 2005) this requirement is relaxed. Finally, there is nothing to prevent the 
scenario described here from creating all the KBOs in situ between 30-44AU, if the solids 
mass were appropriately enhanced and nebula parameters slowly varying (at least for the 
p—1 models), but the cleanup and Neptune migration problems would remain. 

3.5 Other comparisons between model predictions and observa- 
tions 

We have noted several times that the normalized IMFs (figures 4-6) can be made to approach 
or exceed unity (that is, produce the needed mass in planetesimals in the allowed time) if our 
assumed value of tpa = tged, which determines Pgoai is several orders of magnitude too large 
(see for instance Chambers 2010). In this section we discuss other constraints the model can 
be compared with, which are not dependent on this uncertain parameter. 
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3.5.1 Clump encounter times and the age spread in chondrites 

As discussed in section 2.1, recent work suggests a wide spread in the formation times of 
chondrules in a given chondrite - a significant fraction of a Myr. Some concerns remain 
regarding the interpretation of these data as age differences, which will surely be addressed 
as more data emerge. Here we take the age differences at face value, and explore their 
implications for our models. Below we estimate how long a newly formed chondrulc must 
wander through the nebula before, along with many other chondrules being independently 
formed and wandering about, entering one of the rare, dense clumps that is destined to 
become a sandpile planetesimal under the scenario presented here. We would expect this 
timescale to approximate the half- width in the formation age range observed in a particular 
chondrite, by the statistical nature of the process. That is, a few new chondrules will accrete 
shortly after their formation, most will accrete after a time tgnci and some unlucky ones might 
need to wait another tenc or so to find their parent proto-sandpile clump. We assume a fixed 
particle density spatial distribution of proto-sandpile clumps of size I and mass density $ 
even if the individual clumps are everchanging. A wandering chondrule, being nearly tied to 
the turbulent gas, sweeps through the frame defined by the orbiting gas (and dense clumps) 
at a speed Vp roughly equal to the turbulent velocity Vl = coi^l'^ (Cuzzi and Hogan 2003). 
However, because the motion of preferentially concentrated particles is not random in space 
and does not sample all fluid volumes with equal probability, a simple random- walk, particle- 
in-a-box encounter calculation is inappropriate. 

Instead, we use a "duty-cycle" approach similar to that described by CHPDOl (their 
section 6.2). In the Appendix we describe some of the details involved in translating the 
nomenclature of CHPDOl to that used here. CHPDOl integrate a two-dimensional function 
such as our P($,S') over S and distinguish Fv($), the fraction of volume lying in zones of 
mass loading from -Fp($), the fraction of particles lying in such zones. These are not equal 
because particles preferentially are found in dense zones, not randomly in space. CHPDOl 
demonstrated that the cumulative fraction of particles Fp{> $) lying in regions of density 
larger than $ is the same as the fraction of time spent by a given particle in regions with 
density larger than $, Fi{> $). This would be true for the differential functions Fp{^) and 
Fj($) as well. We generalize here to the two-dimensional function Fp($, S) because we are 
more carefully treating the role of enstrophy, but the same identification will hold between 
Fp{^,S) and Ft{^,S), and for their cumulatives. We also adopt a different treatment of 
the cumulative of Fp, as described in the Appendix: specifically, we calculate the fraction of 
particles Fp{> T) lying in zones having properties anywhere within the stable region defined 
by the thresholds Sjnin,^2-, and $i(S') (see figure 3; also see Appendix for derivation of 
Fp{> T)). As in CHPDOl we set Fp{> T) equal to the fraction of time Ft{> T) spent by any 
given particle in zones capable of becoming sandpile planetesimals. We note that Fp(> T) 
is calculated at the level N defining the mode, or maximum, in the IMF for each parameter 
case (figures 3-6), and is thus associated with a lengthscale /. 

For a wandering particle CHPDOl define tin = l/Vp as the time it spends traversing a 
clump of size /, and tenc as the time between encounters with such a clump. Setting tj„ = 
neglects the "peloton effect" in which some (but not most) particles are seen to follow a given 
clump for an extended period of time. Then the duty cycle, or fraction of time spent by a 
particle in regions of size / capable of becoming sandpile planetesimals can be approximated 
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Table 2: Summary of all predictions for the models of figure 4a (relevant to 2.5 AU, and assuming the $i 
of equation 6), as designated by ao,Pgo, and fio- N*, $*, and P* characterize the mode, or peak, of the 
IMF for each case. The encounter time of a chondrule with its ultimate planetesimal-forming clump is also 
tabulated {tenc, in Myr; equation 17), as well as the corresponding value of the normalized radial diffusion 
width Aa/a in a time tenc (sect. 3.5.2), and the implied mass accretion rate M (Mq/jt) given the other 
parameters (sect 3.5.3). We assumed O = 5 x 10~^, appropriate at 2.5 AU. The value of Pgoai assumes 

tpa ~ t'sed' 



(assuming tin < tenc) by 



''enc ''p''enc Lf'enc 



We recall that only the subset of such regions which lie within a fraction of the nebula's 
vertical extent are candidates to become sandpilc planetesimals (sections 3.1, 3.3.1). The 
time fraction spent by particles in this subset of regions is therefore F/(> T) = /5^/^Ft(> T). 
We then obtain the encounter time of a particle with proto-sandpiles by setting tin /t enc = 
Fl{> T) = ^V2i7;(> T) = /3V2Fj,(> T), and solving for tenc- 

I _ 2-^/3gaV^ 2-^/=^ 

- ^^nFp{> T)Vp ~ P^/^Fp{> T)ca^/^ ^ l3^'^Fp{> T)Q' ^ ^ 

In equation (17) above we have used / = 2~^/^L = 2^^^^Ha^/'^ and Vl = ca^^"^, and 
Fp[> T) is evaluated at the value of N giving the peak of the IMF. In the Appendix 
we note that Fp{> T) can be simply related to the modal peak value P*. Of course, these 
arguments are simplified and need to be explored in more detail statistically and numerically. 
Nevertheless, the crude estimates shown for tenc in Tables 2-5, ranging to values of a fraction 
of a Myr, confirm that, in this scenario, accretion is a drawn-out process with timescales 
compatible with those observed. Recall that the observed age dispersion half-widths in 
several different chondrite classes are a few x 10^ years (section 2.1). Cuzzi et al (2010) show 
that the distribution of chondrule ages in two primitive chondrites is compatible with a 
Poisson arrival time distribution characterized by teUC — 0.2-0.4Myr. 



3.5.2 Radial diffusion and "zoning" in the asteroid belt 

Tables 1-4 also show the extent of radial diffusion Aa over the timescale tenci ^ normalized 
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Table 3: See table 2 caption; this table replaces $i of equation (6) by ^sek (section 3.2.1), and assumes 
A/Ao = 10 (sec figure 4b). In this set of models the encounter times are mostly too long, and the radial 
diffusion probably too large, to match the meteoritic constraints. Also, as shown in figure 4b, the desired con- 
dition P* /Pgoal is not generally satisfied. Smaller values oitpa/tsed would increase P* /Pgoal > 1 accordingly 
(section 3.3.1). 
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Table 4: See table 2 caption; this table replaces $i of equation (6) by ^sek (section 3.2.1) and assumes 

A/Ao = 30 (see figure 4c). Here, the meteoritic constraints of age variance (tenc < 1 Myr)and radial gradients 
(Aa/a < 1) may be satisfied by several parameter sets (rows 2-6); note how dramatically the results vary 
from those in Table 3 for only a factor of three change in solids abundance. Figure 4c also shows how the 
P*/ Pgoal > 1 condition is now robustly satisfied. 
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2 4 6 8 

radial distance (AU) 

Figure 7: Diffusion profiles for narrow annular sources at 2AU and 4AU, after 10^, 2 x 10^, 4 x 10^, and 
10^ years. Note that after 10^ years, nearly complete separation is retained between sources at the inner and 
outer edges of the current asteroid belt, but mixing rapidly increases with time and by 10® years the region 
is fairly well mixed. 

by semimajor axis a=2.5AU; a constant ratio H/a — 0.05 is assumed here. The outcome 
of such a random walk in cyhndrical geometry is a shghtly non-gaussian profile centered 
on the starting position, having a halfwidth Aa ~ l.l{Vtencfl'^ where V ~ acH is the 
diffusion coefficient (Cuzzi et al 2003). Figure 7 shows a slightly more detailed treatment, 
in which the constant viscosity Green's function for cylindrical geometry (Cuzzi et al 2003) 
is used to illustrate the diffusive spread of initial delta-functions of tracer "chondrules" 
released at 2AU and 4 AU in a nebula with a — 10~^, after periods of 10^ — 10^ years. 
As indicated by the overlap of the curves in figure 7, some level of discrimination can be 
preserved over timescales of a fewxlO^ years, but mixing is nearly complete by IMyr for 
a = 10~'^. Larger a, of course, leads to more complete radial mixing, as shown in tables 
2-5. This sensitivity of the amount of radial mixing to time, across the 0.1-lMyr range, 
makes emerging developments in chondrule age dating (both the observations themselves 
and the interpretation of the results) highly relevant (section 2.1). 

3.5.3 Nebula mass accretion rate 

Tables 2-5 also show the nebula gas mass accretion rate that would be implied by the various 
adopted parameters: M = STcauT = STcaacH (Lin and Papaloizou 1985). We assume c = 10^ 
cm/s and H{ao) = ao/20 = 2 x 10^^ cm at ao=2.5 AU. The larger values of a, combined 
with the large values of pgo suggested by models where P* / Pgoai > 1, produce mass accretion 
rates substantially larger than regarded as typical for Myr-old protoplanetary disks; a more 
canonical value is a few xlO~^MQ/yr. However, there is an order of magnitude scatter in 
these mass accretion rates which is apparently real (Calvct et al 2000, Hartmann 2005). It 
is interesting that the parameter range giving the most reasonable M for Myr-old disks also 
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satisfies the tenc and diffusion lengtli criteria tlie best (rows 2-3 for A/Ao = 10 and rows 
2-5 for A/Ao = 30). Several combinations of large a and low Pgo can approach or at least 
suggest nominal values of M in the 10~'^ Mq/jic range; however, the values of tenc are very 
long for these, because of the low probabilities of the appropriate clumps, given by Fp{> T), 
and combined with the large a values, radial mixing over these timcscales precludes any 
distinction in properties between contemporaneously formed planctcsimals a few AU apart. 
This constraint only applies to asteroids; no constraint of this type is yet known for KBOs 
where a might well be large. 

3.6 Optimal parameter range 

Looking at all the predictions together {eg. Table 5), one tends to favor the lower values of a 
combined with relatively large gas densities because the combination leads to plausible accre- 
tion rates, Myr-or-less encounter times, and small radial diffusion lengthscales. Somewhat 
larger f3 values, closer to canonical, might be allowable if we adopted a shorter timescale 
for tpa than tged, which brings the normalized IMFs closer to unity by decreasing Pgoai- Ex- 
tremely large radial diffusion ranges tend to characterize the larger values of a ~ 10~^, and 
unless the currently inferred range of chondrule formation ages is misleading due to par- 
ent body resetting (section 2.1), large values of a at 2.5AU are probably inconsistent with 
evidence that ordinary chondrites and CO chondrites - having very different chemical and 
isotopic properties - arc about the same age (Kunihiro et al 2004; Kurahashi et al 2008). 
We note that three grouped, H-like chondrites have been found which are unusual for ordi- 
nary chondrites in having abundant CAIs and CO-like matrix (Kimura et al 2002). These 
chimeric objects may be a sample of a parent body that accumulated at an intermediate 
location between the H chondrite parent (s) and the CO chondrite parent (s), containing a 
blend of components which dominated in the two locations. 

The region of parameter space which is consistent with all the constraints we have men- 
tioned seems to be fairly small {a ~ 10~^, Pgo ~ 1 — 3x10"^, /3o ~ 10~^, and A/A^ ~ 10 — 30); 
in this sense the agreement of the model with expectations is sensitive to small changes in 
model parameters. However, we are trying to match a number of independent observa- 
tions at once, with a very simple model, so we are encouraged that there is any reasonable 
combination of parameters that comes close to matching them all. 

3.7 Caveats and future work 

The models presented here represent the most obvious and straightforward implications for 
primary accretion of the physics of turbulent concentration of small particles, and survival of 
dense clumps of them, as outhned qualitatively in a series of past papers (CHPDOl, CW06, 
CHS08). In order to achieve even this first sanity check, a number of simplifying assumptions 
were made, and a broad range of nebula parameters was sampled. The caveats we feel are 
most important to mention (in rough priority order) are: 

1) Perhaps the most significant "known unknown" is the timescale tpa used to constrain 
the primary mass accretion rate and assess which parameter sets can create the needed 
pre-depletion mass of planetesimals (section 3.3.1). In this paper we assumed tpa = tsed ~ 
100-1000 orbit periods, while Chambers (2010) has assumed tpa — te{l) ~ 0.01-0.1 orbit 
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Table 5: Summary of Model results at 2.5AU, combining A = IQAo and A = ZQAo, from tables 3 and 4. 

periods. An interesting compromise might be tpa ^ ti, or about an orbit period; this might 
be the timescale on which the density field is independently refreshed. A considerable amount 
of study, using 3D numerical models and following clumps for long amounts of time as they 
develop and dissipate, is needed to understand this timescale. 

2) The survival conditions for dense clumps based on our various thresholds (section 3.2) 
he well down a steeply falling slope of the PDFs predicted by our cascade models; these 
cascade models must be checked and tested regarding some of the assumptions built into 
them. For instance, it has been assumed that the multipliers in these models, determined 
near the dissipation scale in numerical models, actually apply over a much wider range of 
scales up the inertial range. While our own tests to date have supported this assumption, 
some other results suggest a scale-dependence to the process. The form and scale-dependence 
(if any) of the multipliers is fundamental (the shapes of the PDFs on their steep edges are 
sensitive to details) and must be checked more closely. The details of the cascade PDFs will 
have implications for other concerns below. 

3) Some of the physical assumptions made by the model (the density of a clump is so high 
that it moves at near-Keplerian velocity) are not fully compatible with the typical values of 
$ ~ 10 that emerge from requiring the IMFs to match the modal asteroid masses and mass 
production rate "goals" (Tables 1-4, figure 4). For $ ~ 10, tged is closer to 1000 orbits than 
100. Fortunately, it seems that a slower contraction poses no obvious problems as long as 
the clump is stable against ram pressure disruption; once it is stable by the criteria of section 
3.1, it only becomes more stable as it shrinks (CHS08). In fact, lower-$ clumps might incur 
lower (3 if their orbital speeds are subkeplerian. However, turbulent eddy variations on this 
long a timescale might assume a larger role in clump disruptions. Of course, higher-^ clumps 
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are also forming, but at lower volume fractions ( eg. , figure 3) . Future studies should explore 
this aspect of the apparent preferred parameter range more carefully with numerical models. 

4) The survival threshold $i (or ^sek), involving centrifugal balance, apphes the inverse 
of an eddy lifetime as if it were a true rotational frequency. This might be overly conservative, 
and if so, the diagonal lines associated with $i (figure 3) shift to the right, increasing the 
values of P($, S) in all the primary IMFs. 

5) The survival threshold $2 involves a poorly-determined parameter (the so-called grav- 
itational Weber number We}.) that quantitatively affects $2 and thus the magnitude of 
P($, S) in all the primary IMFs; numerical simulations of clump disruption must be carried 
out at higher resolution to better constrain Wcq. 

6) Implicit in the primary IMF rates is that all solids at the modeled place and time 
have sizes and densities which are suitable for turbulent concentration; since this is unlikely 
to be the case, some other inefficiency factor must be allowed for, which will increase the 
P($, S) needed to match mass production estimates (which are, of course, quite uncertain 
themselves). As seen in figure 4, order-of-unity changes in A/Ao can provide IMFs where 
P* / Pgoai increases significantly, allowing room for such inefficiency. Moreover, the value of 
tpa/tsed also affects P* / Pgoai (section 3.3.1); our choice of tpa = tged might be relaxed and 
increase P* /Pgoai by 1-3 orders of magnitude, ample to compensate for inefficiency due to an 
initially broad particle size distribution. The recent work of Zsom et al (2010) incorporates 
extensive new experimental results and finds that, under a range of nebula conditions, grain 
aggregates in the asteroid belt region reach a "bouncing barrier" at masses not too different 
from those of chondrule precursors. Such a moderately narrow size distribution of chondrule 
precursors would influence the narrow observed size distribution seen in chondrites, with or 
without turbulent concentration, and fimit the degree of inefficiency inherent in a potentially 
broad pre-TC size distribution. Meteorite data show a range in mean chondrule size across 
chondrite types, even while the shape of the distribution remains apparently invariant (sec- 
tion 2.1). More data of this type is needed for more chondrule types, and for more types 
of objects within chondrites: metal particles, CAIs, and so on, to see how influential the 
limitation of chondrule sizes by limiting the size distribution of their precursors might be. 

7) If Po is actually much lower than the canonical 10~^ in regions where planetesimal 
precursor clumps form (as it appears), the role of turbulent nebula pressure fluctuations as 
an independent disruption mechanism should be reassessed using detailed numerical models 
(section 3.2). Moreover, if this were a global value instead of a being locally-determined hy 
mass loading, as we have suggested, the supply of material from the outer nebula to the 
inner nebula (and thus A/Ao) will be affected (Cuzzi and Zahnle 2004, Ciesla and Cuzzi 
2006). On the other hand, relaxation of our assumption that tpa = tged (section 3.3.1) might 
allow larger values of /3 to satisfy the mass creation rate constraints. We note that recent 
studies of nebulae with "dead" or at least "dull" zones embedded within MRI-active layers 
show higher gas densities than canonical for the terrestrial planet region, and also suggest 
that radial pressure gradients might be more complex than expected from simple powerlaw 
radial dependence (Zhu et al 2010). 

8) We have taken at face value several recent reports of nearly Myr-variance in the ages 
of chondrules found in the same chondrites. If uncontaminated by parent body resetting, 
these variances are critical constraints on primary accretion; a subset of our models indeed 
comes close to explaining this large variance while explaining other chondrite properties. 
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but the parameters are not necessarily in the canonical range. It is important to extend 
relative age measurements on individual objects in the same meteorite to a greater variety 
of samples, and to address concerns that the inferred age differences merely represent parent 
body resetting events (section 2.1) . 

4 Conclusions 

Traditional models of incremental accretion lead to powerlaw size distributions with equal 
mass per decade for particles ranging between centimeters and (at least) tens of km in size, 
where runaway gravitational growth sets in; such distributions are increasingly thought not 
to be compatible with the current asteroid size distribution (Bottke et al 2005, Morbidelli 
et al 2009a, Weidenschilling 2009). We have used a very simplified physical model to pursue 
the most obvious implications of the fate of dense clumps of mm-sizc particles, which arc 
aerodynamically selected for preferred concentration in nebula turbulence. We follow the 
simplest physics determining which dense clumps avoid disruption and evolve, on periods 
of 100-1000 orbits, into objects with some physical cohesion. We find that some small 
fraction of these dense clumps (those forming near the nebula midplanc) can proceed to 
become "sandpile" planctcsimals having diameters in the 20-200 km range observed for 
today's asteroids, in an abundance which is broadly consistent with poorly known estimates 
of the mass of the pre-depletion primordial belt. That is, it appears possible for most of the 
mass in primitive bodies to have simply skipped over the problematic m-km size range. 

The model is a simple one; except for the physics and statistics of the turbulent cascade 
model, it contains little more than physical scaling arguments. The parameter range studied 
is broad but not exhaustive, and the preferred range is not optimally centered on "canonical" 
nebula conditions. Specifically, lower-than-expected headwind magnitudes, and higher-than- 
canonical gas densities and solid/gas ratios, are needed to achieve quantitative matches to 
(uncertain) estimates of required mass production rates (section 3.4). We have provided 
some thoughts on why these conditions might not be unrealistic in view of the uncertainties, 
but more work along these lines is surely needed. For instance, the critical timescale tpa 
remains poorly understood even in principle (sections 3.3.1 and 3.7; Chambers 2010). 

In spite of the fact that the optimum parameter set is perhaps not the canonical one, 
we are encouraged that such a simple model can achieve so much with any plausible range 
of parameters. The model scenario potentially explains many things simultaneously: the 
mean size, the size distribution, and the age dispersion of chondrules in chondrites, and the 
modal mass and mass production rate of primary asteroids. While doing this, the scenario 
also provides a plausible radial mixing length for planetesimal constituents between their 
creation and their accretion, allowing for some radial zoning of primitive asteroids, while 
retaining a canonical nebula mass accretion rate for Myr-old protoplanctry nebulae (for 
which considerable scatter does exist, however). This is not to say that we believe a final 
solution is at hand; there are many serious uncertainties and unresolved issues (section 3.7). 
Moreover, this scenario may not work alone; planetesimal formation is likely to have been 
very complicated. Aggregation effects acting on small particle scales before chondrules were 
even melted (most recently Giittler et al 2010 and Zsom et al 2010) and particle-gas dynamics 
acting on larger particles than we consider (Johansen et al 2007) both operate in the same 
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(turbulent) environment we have studied. 

The scenario presented here identifies a path which - while still fraught with hazards - 
leads directly from freely-floating, mm-size nebula particulates, to sizeable (tens to hundreds 
of km diameter) sandpilc planetesimals formed almost entirely of size-sorted particles, which 
is in reasonable accord with the meteorite record. The typical encounter time of any par- 
ticular chondrule with the clump in which it becomes a planetesimal can be a significant 
fraction of a Myr - comparable to formation age variance observed isotopically in several 
different chondrite groups. The constituents can be narrowly sorted physically because pri- 
mary accretion refiects local properties, but can be diverse chemically and isotopically, with 
well-defined group properties only in the ensemble, because production/alteration regions 
can be separated in space by several H, and in time by nearly 1 Myr. The primary sand- 
piles, made of constituents with a wide range of formation ages and chemical and isotopic 
properties, are in the size range of today's "fossil" asteroids, and we imagine they go on to 
experience an extended sequence of compaction, heating, and sintering, perhaps even before 
the dynamical depletion stage when more violent, erosive and destructive impacts lead to 
the objects we see today. The process is capable of starting very early in nebula history and 
proceeding for a long time, as nebula solids evolve. Primary objects of tens-hundreds of km 
diameter forming early will almost certainly melt extensively, while those forming more than 
1.5-2.5 Myr after CAIs will be able to remain unmelted and "primitive". 

In this scenario, one would expect some degree of both radial mixing and temporal evolu- 
tion in the properties of chondrulcs ending up in a particular chondrite. Along these lines it 
is useful to recall that, while the properties of a chondrite group are well-defined, there is not 
only the formation age diversity mentioned above, but also a substantial variance in petro- 
logical, chemical, redox, and isotopic properties amongst the constituent chondrules in any 
single chondrite (Scott and Krot 2005; Brearlcy and Jones 1998). The well-defined properties 
of a chondrite group may only manifest the ensemble homogeneity of their common parent 
body - which itself might represent a grab-sample of constituents which came together at 
a given time and place. Other planetesimals forming nearby and contemporaneously might 
have very similar properties, but planetesimals forming at a different time, or in a differ- 
ent location, would draw from a slightly evolved mixture of essentially the same building 
blocks, perhaps modified by ongoing alteration, remelting, changing oxidation environment, 
and/or mixture with particles from adjacent regions, resulting in different ensemble proper- 
ties. Cuzzi et al (2005) discuss other meteoritics implications in more detail, including the 
concept of "complementarity" . Some models of the subsequent stage of accretion, charac- 
terized by gravitational scatterings and collisional mergers of primary objects, suggests even 
more radial diffusion (Bottke et al 2006, Levison et al 2009). 

As a cautionary remark, we note that the meteorite data, interpreted in the context of 
our scenario, also suggest that combinations of nebula pg and a must have varied in time 
and/or space over the region and duration of primary accretion. This is implied because 
different chondrite groups with similar accretion ages (ordinary and CO chondrites) have 
noticeably different modal chondrule sizes - varying by a factor of several (King and King 
1978, Rubin 1989, Scott and Krot 2005; Brearley and Jones 1998; Kurahashi et al 2008). In 
the context of turbulent concentration, because all chondrules are after all made mainly of 
silicates and have comparable densities, variations in a and pg are the most obvious way to 
do this (CHPDOl, figure 1). 
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Extension of the scenario to the Kuiper Belt region has also been explored; IMFs show 
a similar preference for 10-lOOkm diameter objects. Mass creation rates assuming canonical 
cosmic abundance tend to fall short of estimates for KBOs (which are rather poorly known, 
however), but enhanced solid abundances over cosmic dramatically improve the agreement 
with mass production estimates. Moderately flat nebula gas density distributions are sub- 
stantially more favorable to the extension of this scenario from the asteroid belt to 30AU. 
Under such conditions, and given the strong dependence of primary accretion on local solids 
abundance, the diversity of KBOs - from thoroughly melted, water-ice-mantled objects such 
as Haumea to those retaining abundant "supervolatiles" - might be explained by the same 
drawn-out accretionary process as we envision for the asteroids. The specific scenario ex- 
plored here assumed KBO initial formation between 16-30 AU and subsequent dynamical 
emplacement to 30-44AU; this is not a requirement of our model, however. Crude scahng 
estimates suggest that the mass of the "cold classical" KBOs could be formed in situ, by the 
physics discussed here, from a local mass density that would not lead to excessive migration 
of Neptune or a problematic subsequent cleanup. More refined future models of planetesi- 
mal formation in the 16-30AM region should, for self- consistency, strive to produce roughly 
twice the mass we assumed for this first assessment, to allow for the cores of the ice giants 
themselves. 
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Appendix 

Here we provide some notational clarifications to better connect the distribution func- 
tions P($, S) and Fp{^, S) of this paper to similar functions in CHPDOl and Hogan and 
Cuzzi(2007). We first note that we have defined P($, S) in this paper as a probability per 
unit logio($) fT-iifl logio('S') (as in Chambers 2010; a hand check for identical parameters 
shows that our PDF contours are in very good agreement with those of Chambers). Hogan 
and Cuzzi (2007) are not specific about the use of log^Q vs. logg = In, but here we are 
more explicit. We define P'{^,S), which Hogan and Cuzzi (2007) call P{^,S), as a true 
differential probability density per unit $ and per unit S. Since, for instance, S — e^°'^, 
logS — loge • \nS and dlogS = loge ■ dS/S. Wc require the functions P($, S) and P'($, S) 
to be separately normafized: /~ P'($, S)d^dS = 1 = /o°° P{^, S)dlog^dlogS . Since 
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P'($, S)d<^>dS = P'($, • dlog<^/\oge){S ■ rflog5/loge), it follows that 

The quantity P'($, S)^S was proposed by Hogan and Cuzzi (2007) as a convenient estimate 
of effective volume fraction "at" S), as it covers the range ($ ± A$/2, S ± AS/2) with 
A$ = $ and AS* = S, but here we use instead the quantity P($, S) in this role. The 
difference is a factor of log^e, and can be explained by P($, 5*) having a larger effective 
binning size A$ = ^/loge, AS = S/loge. In fact, we find that P($, 5) actually better 
approximates the more formally exact integrals (below) than does P'($, S)^S. 

To be even more specific, we define P($, S) by taking the number of cascade outcomes 
H which lie in bins at ($, S) within some bin size dlog^ ■ dlogS = 5^, and normalizing by the 
total number n of cascade outcomes so P($, S) = H{^, S, d)/nd^. This makes P($, S) also 
a probability (fractional volume) density, but per unit log^ol*^)) P^r unit logiQ(S'). It is easy 
to show that using this definition and P'($, S)^S = log^eP($, S), P($, S) and P'($, S) are 
both automatically normalized: P'($, S)dMS = P($, ^)rflog<l> ■ dhgS = 1. 

We also note that the particle concentration C of CHPDOl is simply related to the mass 
loading factor $ of this paper and CHS08. Specifically, C = Up/ (up) = Pp/ (pp), where Up is 
the number of particles per unit volume, (n^) is its global average, and all the particles are 
of equal mass in our simple treatment. Thus $ = pp/pg = A'C, where A' — 0.01^4/740. As in 
CHPDOl we must distinguish between the fraction of volume found at some ($, S) and the 
fraction of particles found there, because the particles are not randomly distributed in ($. 5") 
space. Because the volume and particle fractions of CHPDOl (Pi/(C), Fp{C)) are normalized, 
they can also be written as Py($),Pp($). Then recaU that Py($) = P'{^,S)dS so 
Fv{> $) = P'{^',S)d^'dS ^ 1 as $ ^ 0, as in CHPDOl. We define a particle 

fraction Pp($, S) oc $P'($, 5"), which is the 2D extension of the function Fp{C) oc CPy(C) 
of CHPDOl, and like it, must be separately normalized. We define a normalization constant 
Cp such that Cp $P'($, S)d^dS = 1; it can easily be shown that the same Cp implies 

that Cp $P(<I>, 5)rflog$ ■ dlogS = 1. 

The closest analog to Fp{C) of CHPDOl is Pp($) = Cp$Pi/($) = Cp$ P'($, S)dS, and 
has the equivalent cumulative function Fp{> $). However, for the purpose of the current 
paper we are working in a more profoundly 2D regime, where both $ and S are important. 
Thus we will work with the cumulative of Pp($, S) = Cp^P'{^, S) that represents the fraction 
of particles lying in clumps that exceed all of our thresholds ($i, $2, 'S'^m) for sandpile 
formation or primary accretion. Any clump lying above the entire Smin — ^1-^2 threshold line 
is capable of becoming a sandpile (section 3.3.1), and encounters of a wandering particle with 
all of them should be included in estimation of ignc (section 3.5.1). Thus we define Pp(> T) - 
the particle fraction lying in all proto-sandpile-clumps at any given time - as an integral over 
the 2D segment of ($,5*) space lying above the threshold value ^t{S) = ■max{^2,^i{S)) 
and to the right of Smin (see figure 3). That is, 

TOO ^00 ^pu^ s)d^dS r°° r°° 

Fp{> T) = 'tJ^^I^,,^ ,^ ,^ = Cp / / ^P\^,S)dMS, (18) 

Jo Jo ^P'{^.S)dMS JSmiuUriS) 

where the constant Cp = Ao/A. This cumulative measure will be dominated by clumps falling 
closest to the threshold line, and in particular at the peak of the IMF (figure 3) which we 
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characterize by $ = $* and P($*,5'*) = P*. For example (see section 3.3.1 and Chambers 
2010) the total primary accretion rate of planetesimals Mp^ can be written as 



where the double integral is over the entire range of jS" where sandpiles can form, the 
integrand is the product of particle mass per unit volume in a clump $pg times the volume 
fraction in clumps^ and the numerator is the volume in which this transpires using only a 
narrow region of thickness /3^/^ near the midplane. The denominator tpa is some formation 
timescale, which we take as tged but could be shorter (section 3.3.1; Chambers 2010). We 
can rewrite this as 



pa J Smin 

From expressions given above, it is clear that , J^^^g^^P'{^, S)d^dS — Fp{> T)/cp — 
Fp{> T)A/Ao. Meanwhile we could also approximate Mpa ^ 

^pa 

where the peak of the IMF is at (A^*, $*, S*) and has associated P($, S) = P*. Physically this 
amounts to saying that the bulk of the primary accretion is that which occurs near the peak 
of the IMF and neglecting the contributions from further down the IMF, but treating the 
volume fraction as a binned value over a range A$* = $*/loge, AS* — S*/loge as described 
above, based on our definition of P{^,S). Comparison of the above equations shows that 
taking $*P*($*,S'*) = Pp(> T)A/Ao makes the expressions equal. Wc make use of this in 
section 3.3.1, using numerical vahdation (in tables 1-4) that $*P*($*,5'*) = Fp{> T)A/Ao 
is indeed valid to tens of percent. As described in section 3.5.1 and CHPDOl, we can also 
associate Pp(> T) with the fractional time Pt(> T) spent by a given particle in protosandpile 
clumps. 
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